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Abstract 



' We propose and analyze numerical methods for the Heath-Jarrow-Morton (HJM) 

■ model. To construct the methods, we first discretize the infinite dimensional HJM 
equation in maturity time variable using quadrature rules for approximating the 
arbitrage- free drift. This results in a finite dimensional system of stochastic differ- 
ential equations (SDEs) which we approximate in the weak and mean-square sense 
using the general theory of numerical integration of SDEs. The proposed numer- 
ical algorithms are highly computationally efficient due to the use of high-order 

■ quadrature rules which allow us to take relatively large discretization steps in the 
maturity time without affecting overall accuracy of the algorithms. Convergence 

!y-^ | theorems for the methods are proved. Results of some numerical experiments with 

C\| 1 European-type interest rate derivatives are presented. 

Q\ • Keywords. Infinite dimensional stochastic equations, HJM model, weak approx- 

imation, Monte Carlo technique, interest rate derivatives, method of lines, mean- 
square convergence. 
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x. 

c3 1 Introduction 

The framework proposed by Heath, Jarrow and Morton [13] - HJM henceforth - models 
the evolution of the term structure of interest rates through the dynamics of the forward 
rate curve. These dynamics are described by a multifactor infinite dimensional stochastic 
equation with the entire forward rate curve as state variable. Under no-arbitrage condi- 
tions, the HJM model is fully characterized by specifying forward rate volatility functions 
and the initial forward curve. The original HJM framework is used for modelling fixed in- 
come markets (see [13j HI [6j [7] and also references therein). Recently, the HJM philosophy 
has been extended to credit and equity markets (see, e.g. the recent review [5]). 

The HJM model has closed-form solutions only for some special cases of volatility, and 
valuations under the HJM framework usually require a numerical approximation. As far 
as we know, the literature on numerics for the HJM model is rather sparse. The common 
approach (see, e.g. [12j [15j EJ [3] and the references therein) is to take coinciding grids 
in the running time t and in the maturity time T. The known methods differ in the way 
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they approximate the integral in the arbitrage-free drift of the HJM model while they 
all use Euler-type schemes for discretization in time. In [T2| [T5| [9] approximations of 
the arbitrage-free drift are chosen so that the overall discrete approximations of the HJM 
equation preserve the martingale property for the discretized discounted bond process. A 
different numerical approach, based on a functional backward Kolmogorov equation, is 
considered in [17] . In comparison with other works, the papers (3j [17] rigorously prove 
weak convergence of the proposed numerical methods. 

In this paper we propose and analyze a new class of effective numerical methods for 
the HJM equation exploiting the idea of the method of lines. These methods can be used 
for simulating HJM models of various specifications. Our main focus is on the weak- 
sense numerical methods which can be used for valuing a broad range of interest rate 
products. To construct the numerical methods, we first discretize the infinite dimensional 
HJM equation in maturity time variable T using quadrature rules for approximating the 
arbitrage- free drift. This results in a finite dimensional system of stochastic differential 
equations (SDEs). As we show in the paper, if we take a quadrature rule of order p, the 
solution of this finite dimensional system of SDEs converges to the HJM solution with 
mean-square order p in the maturity time discretization step A. From the method of lines 
point of view, we interpret the maturity time T as a "space" variable and the running 
time t as a "time" variable. To get fully discrete methods (discrete in both T and t), 
we approximate the obtained finite dimensional system of SDEs in the weak and mean- 
square senses using the general theory of numerical integration of SDEs (see, e.g. [18J). 
The proposed numerical algorithms are computationally highly efficient due to the use of 
high-order quadrature rules which allow us to take relatively large discretization steps in 
the maturity time without affecting overall accuracy of the algorithms, i.e., the number 
of forward rates we need to approximate at each time moment t is significantly less than 
it is usually required in the case when the time grids for t and T coincide. Further, since 
we exploit the method of lines, we have flexibility in choosing appropriate approximations 
in "space" and "time" separately. As we will see, in practice it is beneficial to use higher 
order rules for integration with respect to maturity time T and lower order numerical 
schemes for integration with respect to time t. 

The rest of the paper is organized as follows. In Section [2] we recall the HJM frame- 
work. Section |3] deals with construction of a class of numerical methods of a general 
form for the HJM equation. We start with T-discretization of the HJM equation (Sec- 
tion [3Jjj, then consider discretization in time t (Section I3.2[) . and finally, based on the 
results of the previous subsections, we obtain approximations suitable for evaluating prices 
of European-type interest rate contracts under the HJM framework (Section 13 .3j) . In Sec- 
tion H] we prove convergence theorems for the methods constructed in Section [3J We first 
prove convergence theorems for the HJM approximation discrete in the maturity time T 
only (Section 14. ip . Then we analyze weak convergence of fully discrete methods to the ap- 
proximations discrete in the maturity time (Section 14. 2p . We show that this convergence 
is uniform in the maturity time discretization step A in order to obtain weak convergence 
of the fully-discrete numerical methods to the solution of the HJM equation. We note that 
both the considered class of numerical methods and proof of their convergence include 
the known numerical schemes for HJM as, e.g. those from [9]. In Section [5] we illustrate 
the introduced class of numerical methods from Section [3] by presenting some particular 
algorithms of various accuracy orders, which are ready for implementation. In Section [6] 
we propose a fully-discrete mean-square approximation of the HJM equation and prove 
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the corresponding mean-square convergence theorem. In Section [7] we test the proposed 
numerical algorithms on the Vasicek and proportional volatility models. 



2 HJM framework of the instantaneous forward rate 
dynamics 

Throughout this paper we assume that there exists an arbitrage-free market with a fric- 
tionlessly traded continuum of default-free zero-coupon bonds {P(t,T), t < T, T £ 
[to,? 1 *] , t £ [i ,£*]}, where P(t,T) denotes the price at time t of a bond with matu- 
rity T. We require that P(T,T) = 1 and P(t,T) is sufficiently smooth in the maturity 
variable T. 

A convenient, albeit a theoretical concept, the forward rate f(t, T),t <T,T £ [to, T*] , 
t £ [to)**]) represents the instantaneous continuously compounded interest rate prevailing 
at time t for riskless borrowing or lending over the infinitesimal time interval [T, T + dT] . 
The relation between zero-coupon bonds and instantaneous forward rates is given by 



P(t,T)=ex P (^ T Jf(t,u)d u y (2.1) 



The current instantaneous rate, or so-called short rate, is 

r(t) := f(t,t). (2.2) 
To represent the accumulating factor one can define the savings account 

B(t) = exp r(s)d^j . (2.3) 

The HJM framework [13] models the dynamics of the forward curve 
{f(t,T), t<T, T£ [t ,T*], te[t ,t*}}. 
Given an integrable deterministic initial forward curve 

f(t ,T) = f (T), 

the arbitrage-free dynamics of the forward curve under the risk-neutral measure Q asso- 
ciated with the numeraire B(t) are modelled through an Ito process of the form 

f(t,T)-f (T) = j\ T (s,T)^j\(s,u)du^ds (2.4) 

t 

+ 1 a J (s,T)dW(s), t <t<t*AT, t <T<T*, 



[ a T ( Sl T)dW(s 

Jt 



where W{t) = (Wi(t), . . . , Wd(t)) is a d- dimensional standard Wiener process defined on 
a filtered probability space J 7 , {^ r i} to<i<i » , Q) satisfying the usual hypothesis; cr(t,T) 
is an Revalued JVprogressively measurable stochastic process with f£ \a(s, T)| 2 ds < oo; 
and t* AT := min (t*,T). 
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In general, the volatility o~(t,T) := <j(t,T,u) can depend on the current and past 
values of forward rates. In this paper we restrict ourselves to the case in which a depends 
on the current forward rate only, i.e., 

a(t, T) := {a x {t, T, f(t, T)), . . . , a d (t, T, f(t, T))) T , (2.5) 

where <7j(t, T, z), i — 1, . . . , d, are deterministic functions defined on [t , t*} x [t , T*] x R. 
Then the term a(s,u)du in (12.41) can be written as a(s,u, f(s,u))du, and, conse- 
quently, (12. 41) -( 1231) is an infinite-dimensional SDE. We impose the following assumptions 
on the HJM model (12^) - (1231). 

Assumption 2.1 The functions <Jj(t, T, z), i — 1, . . . , d, are uniformly bounded, i.e., there 
is a constant C > such that 

\<Ti(t,T,z)\<C, (t,T,z) e [t ,t*] x [t ,T*]xR. (2.6) 

Assumption 2.2 For sufficiently large pi, P2 > 1, the partial derivatives 

dPdTW >Q<J<P^Q<k + l<P2, i = l,...,d, (2.7) 
are continuous and uniformly bounded in [to,t*\ x [to,? 1 *] x M. 

Assumption 2.3 The initial forward curve fo(T), T E [t ,T*] , is deterministic and 
sufficiently smooth. 

The imposed conditions are sufficient to ensure that the SDE ( I2.4l) -( f2~31) has a unique 
strong solution f(t,T), which is sufficiently smooth in the last argument (see [131 EI] an d 
also [HI [8] for differentiating SDE solutions with respect to a parameter). Further, it is 
not difficult to show that they imply boundedness of exponential moments of f(t, T), i.e., 
for a c G M there is a constant C > such that 

Eexp{c\f{t,T)\)<C (2.8) 

for all (t, T) G [to,t*] x [to, T*] . The constant C in (12. 8p depends on the initial forward 
curve foiT), volatility o~(t,T,z), and on c. 



Remark 2.1 As it was shown in \2J$ , for the SDE l \2.4\) - i\2.5\) to have the unique strong 



solution it suffices to require a weaker assumption than Assumption 2.1: 

Wi(t,T,z)\ <c(l + |z| 1/2 ). 

However, in the paper we restrict ourselves to the stronger set of conditions which allow 
us to consider methods of higher order. Assumptions 2.1-2.3 are sufficient for all the 
statements in this paper. The choice of p\ and P2 depends on a particular algorithm ( as 
usual, the more accurate an algorithm the more derivatives are needed). At the same time, 
the imposed conditions are not necessary and the proposed numerical methods themselves 
can be used under broader assumptions. 
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We pay attention that Assumptions 2.1-2.3 do not guarantee positiveness of f(t,T) 
which could be a desirable property taking into account the financial context of the HJM 
model. One can notice that if we also require that 

/ (T) > and a^t, T, 0) = 0, i = 1, . . . , d, (t, T) e [t , **] x [t , T% 

then the forward rates are nonnegative f(t,T) > for all (t, T) 6 [to ; ^*] x [^cb^*]- 

Our main objective is to propose efficient numerical methods for pricing interest rates 
derivatives. Among these instruments are interest rate caps, floors, and swaptions [U El 
|22| 123} [7]. A cap price is obtained by summing up the prices of the underlying caplets. 
Consider a caplet set at time s k with payment date at Sj > s k , with strike K and unit 
cap nominal value. Its price at time to < s k is given by 



Eexp(^-J^ r(u)dv^j 1 - (1 + K(si - s k )) exp ^- f(s k , u)dv^j 



(2.9) 



Now consider a payer swaption of maturity sj~ and with underlying swap maturity Sj > s k . 
Its price at time to < s k can be found as 

Eexp^—J r{u)dv^j 1 — exp ^— J f(sk,u)dv^j (2.10) 
-K ^ (sj - Sj_i) exp f - / f(s k ,u)duj 

j=k+l ^ Js k /. 

In (12 .9p and (I2.10p and in what follows expectation E (•) without any index means expec- 
tation taken with respect to the risk-neutral measure Q. 

Let G(z), z G M, be a payoff function satisfying the global Lipschitz condition, i.e., 

\G(z) -G(z')\ < K\z-z'\, z,z'eR. (2.11) 

In this paper, motivated by the above examples, we consider the price of a generic interest 
rate contract under risk-neutral measure of the form 

F{t , f (•) ; s k , Si ) = E exp(-Y(s k ))G {P{s k , s t )) , (2.12) 

where 

Y(s k )= / r{u)du, (2.13) 

P(s k ,Si) = exp(-Z(s k ,Si)) , (2.14) 

and 

Z(s k ,Si)= f(s k ,u)du. (2.15) 

We note that ( 12. 12ft does not cover the case of swaptions ( 12. lOj) . To include swaptions, 
the payoff G in ( 12. 12ft should be of the form G (P(s k , s k +i), • • • , P(s k , Sj)) and 

F(t ,f (-);s k ,s k+1 ,...,Si) = E exp(-Y(s k ))G (P(s k , s k+1 ), . . . , P(s k , Si)) . (2.16) 

We limit ourselves in the paper to the payoff of the form (" 12 . 1 2[) for the sake of transparent 
exposition. All the proposed numerical algorithms are applicable to the more general 
form of the payoff ( I2.16p . Also, no additional ideas are required to extend our theoretical 
analysis to the case (I2.16p . 
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Remark 2.2 (Forward measure pricing) The HJM dynamics can be written under the 
Sk-forward measure (see details e.g., in J]^\2M\El) instead of the risk-neutral measure. In 
this case the corresponding SDE has the form (cf. flff.^p).- 



t / pT \ rt 

T , 



f(t,T)-f (T) = / a l (s,T)( a(s,u)du\ds+ a 1 (s,T)dW 8 "(s), (2.17) 

JtQ \Js k J JtQ 

to < t < t* AT A s fc , t <T<T*, 

with W Sk (s) being a d-dimensional standard Wiener process under the Sk-forward measure 
Q Sk . The pricing formula for a generic interest rate contract with payoff G (P(sk, Sj)) 
under Q Sk is (cf. flgQg) ): 

F(t , f (•) ; s k , Sl ) = P(t , s k )E Sk (G (P(s k , Si ))) . (2.18) 

T/izs /orm computationally simpler than since it does not require evaluation of 

the short rate. At the same time we note that pricing of some interest rate products 
(e.g., Eurodollar futures) require the use of risk-neutral measure ^ WT^. In this paper we 
construct numerical algorithms for approximating pQJj). Obviously, these algorithms are 
readily (actually more easily) applicable to < \2.18b . 



3 Numerical method 

In this section we construct a numerical method for simulating ( 12.12}) with the forward 
rates f(s k ,u) satisfying the infinite-dimensional SDE f l2.4l) -( l2~5l) . Examples of some par- 
ticular algorithmic realizations of this method are given in Section [5j 

This section is organized in the following way. We first introduce a maturity time 
discretization (T-discretization) and arrive at a finite dimensional approximation of (12. 4ft - 
(12.51) . i.e., at a finite system of SDEs (Section 13.11) . Then (Section 13. 21) we discretize 
time (t-discretization) and apply a weak-sense numerical integrator to the obtained finite 
system of SDEs. Finally, Section 13.31 deals with approximating the functionals Y and Z 
from (jgJSD - flHTgp and the option price (l212|) . 

For the simplicity of presentation, we consider equally-spaced grids for maturity time 
T and time t. A nonuniform discretization might be needed in practical financial ap- 
plications, and a generalization of the proposed algorithms to nonuniform time grids is 
straightforward. 



3.1 T-discretization 

Consider a uniform partition of the maturity time interval [toj^*] with a maturity time 
step (T-step) A = (T* — t )/N : 

t = T <--- <T N = T*, Ti = iA, i = 0,...,N. (3.1) 

Introduce the index notation which we will use throughout the paper. Denote by £(t) 
the auxiliary index dependent on time t so that 

£{t) = max{z = 0,1,..., iV : (3.2) 
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and by g(t) the auxiliary index dependent on time t so that 

g{t) = mm {i = 0,1,..., N : t < T { } , (3.3) 

i.e., Te(t) <t< and (or T Q ^) is the closest node on the grid ( 13. ip to the time t 
from the left (or from the right). We also note that g(t) = £(t) + 1. 

Further, we require for simplicity that A is sufficiently small so that a number of 
nodes Tj between t* and T* is sufficient for realization of all the quadrature rules and 
interpolation/extrapolation used in the method which we introduce in this Section [3j We 
will pay attention to the required amount of nodes between t* and T* in the method's 
description. At the same time, if in practical realization the distance between t* and T* 
is relatively small in comparison with the chosen T-step A, then one would need to run 
simulation for a slightly longer maturity-time interval, extending it beyond T* by a few 
steps of A (see further explanation in Section 15. 3p . 

For a node T i7 i = 0, . . . , N, on the maturity time grid (13. ip . we approximate the 
integrals in (12. 4ft : 



Ij(s, Ti) := I (Tj(s, u)du, j = 1, . . . , d, t < s < t* A Ti, i = 1, . . . , N, (3.4) 

J s 

by a composite quadrature rule Sj.(s, T i} A) : 

«(s,Ti) 

I j (s,T i )&S Ij (s,T i ,A) = A lk(s)°j(s,T k ), (3.5) 

k=g(s) 

where the quadrature rule's weights 7*.(s) and the nodes k = g(s), . . . , k(s, Ti) are chosen 
so that under Assumptions 2.1-2.3 the approximation is of order 0(A P ) for a given p > 1, 
i.e., the numerical integration error is estimated as 

/ r i2\ 1/2 

f^^.^T^A)-/,-^^)] J <CA? (3.6) 

with a constant C > independent of A, s, Ti, j. Some examples of such quadratures 
are given in Section |5l We note (see details in Section [5]) that when s and T« are close, to 
approximate Ij(s,Ti) with a required accuracy the number k(s,Tj) in ( 13 .5p can be chosen 
larger than i. We recall that since we assumed that there is a sufficient number of nodes be- 
tween t* and T* the number k(s, Tj) does not exceed N. We will also use the vector notation 
I(s,T t ) := (J 1 ( S ,T i ),...,J,( S ,T i )) T and^(s,T 4 ,A) := (S h (s, T t , A), . . . , S Id (s, T h A)) T . 

For a fixed T — T%, it is convenient for later purposes (namely, for computing the short 
rate r(t) = fit, t) as it will become clear in Section [373]) to consider the SDE (I2.4p -( l23]) on 
a slightly larger time interval: t Q < t < t* A T( i+ i) AN , i.e., for Tj < t* (note that t* < T N ) 
we would like to extend the definition of f(t,Ti) from t e [tcb^i] to t € [t Q ,T i+ i\. Though 
from the point of view of financial applications the forward rate f(t,Ti) is not defined 
on the interval t e (Tj,Tj + i], Assumptions 2.1-2.3 guarantee that (I2.4p -( l231) has the 
strong solution on the extended interval and, as it will be seen in future, this extension 
is beneficial from the computational prospective (see also Remarks 13.11 and 13. 2j) . This 
extension requires from us to consider, in addition to (13 .4p . the integrals 

rT l(s) 

Ij{s,T e{s) ) := / a j (s,u)du. (3.7) 
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We approximate these integrals by a quadrature rule analogous to the one in ( 13 .5p but 
with summation index k starting from £(s) : 



I j (8,T tw )*S Ij (s,T t{B) ,A) = A lk{s)a 3 { S) T k ), (3.8) 

k=l(s) 

and we require that its error satisfies (13. 6p . Combining (I3.5P and (I3.8p . we will write in 
what follows that 

K(s,Ti) 

S Ij (s,T i ,A) = A ]T lk (s)a j (s,T k ) (3.9) 

fc=l(s) 

with the coefficient j£( s )(s) = if i > £(s). 

Using (13. 9p . we approximate the solution f(t, T) of the infinite-dimensional SDE (12. 4ft - 
(12 .5p at the nodes T — T , . . . ,T N , by the N + 1-dimensional stochastic process f l (t) « 
f(t,Ti), i = 0, . . . , N, which satisfies the finite system of coupled SDEs: 

r(t)-f = fa] \s)~Sj{s, T i ,A)ds+ [ a](s)dW(s), (3.10) 

J to J to 



to J to 

t < t < t*AT (i+1)AN , i = 0,...,N, 



where 



Po = fo(Ti), (3.11) 
5»(s) = (a i:1 (s) } ... } a i4 (s)) T = (ax(s,Ti,f (s)), . . . ,a d (s,T i7 f (s))) T 

§i(s, T h A) = (s h (s,T i ,A),...,S Id (s,T i ,A)Y 
and 

K( S ,Ti) 

Si^s, Ti, A) = A^ lk (s)a kd (s). (3.12) 
We emphasize again that we extended the time interval from t G [£o ; £*] to t £ [to,t* A 

T(i + l)AN]- 

Assumptions 2.1-2.3 guarantee the existence of the unique strong solution of (I3.10p - 
(I3.12p . Further, it is not difficult to show that they also imply boundedness of exponential 
moments of /*(£), i.e., for a c £ M there is a constant C > such that (cf. (12.81) ): 

Eexp(c\P(t)\)<C (3.13) 

for all t £ [t , t*) A T(j +1 ) AA r, i — 0, . . . , N. In connection with (13.131) we recall that due to 
Assumption 2.3 the initial forward rate curve fo{T), t <T < T*, is bounded by a finite 
constant. Hence, / l (0) = /q, i = 0, . . . , N, are bounded by the same constant. 

In Section HJ] we prove (see Theorem 14.11) mean-square convergence of f l {t) to f(t,Ti) 
when A — > 0. We note that the system (I3.10p - (l3.12p plays only an auxiliary role in our 
consideration. It is used as guidance in constructing fully discrete numerical algorithms 
(i.e., discrete in both T and t) and also in proofs of their convergence. 
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3.2 ^-discretization 



In this section we discretize the finite system of coupled ordinary SDEs ( I3.10l) - (l3.12p with 
respect to time t and thus arrive at a fully discrete method. 

We introduce an equally-spaced grid for time t with step (t-step) h = (t* — to)/M: 

t < ■ ■ ■ < t M = t*, t k = kh, k = 0,...,M. 

In what follows we use the notation (cf. CI3 . 21) and (13.31) ) : 

4 := £(t k ), Q k ■= Q(t k ). (3.14) 

We consider an approximation fl +1 of f l (t k+1 ) from (13.101) (i.e., a full discretization 
of ([23D-(E3]) in both T and t) of the form 

/* = /*, i = 0,...,N, (3.15) 

fk+i = fk + A %, Tf, fi, j = 4+i, • • • , «(**+!, Ti) V i; h; 
i = £ k+ i,...,N, k = 0,...,M, 

where the form of the functions A 1 depends on the coefficients of (13. 10[) -( T3. 12[) . i.e., on 
a and on a choice of the quadrature rule Sj .; n(t k ,Ti) is as in the quadrature (13.121) ; 
£ k , k = 0, . . . , M, are some random vectors which have moments of a sufficiently high 
order and £ fc for k > are independent of /j, i = £j, . . . N, j = 0, . . . , k, and of £ , . . . , Cfe-i- 
To simplify the exposition of our theoretical analysis, in what follows we consider the 
extended f l (t) and f k . We put 

f (t) = f (T i+1 ), T (i+1)AN At*<t<t*, 0<i<£(t*)-l, 

and then the N + 1-dimensional vector {f l {t), i = 0, . . . , N} is defined for all t G [to, £*]• 
We put 

= fc = m + l,...,M, 0<i<*(O-l, 

where m = [(Tj + i —to)/h~\ — 1 (we recall that [•] denotes the ceiling function). Then 
the iV + 1-dimensional vector {fl, i = 0, . . . , A^} is defined for all k = 0, . . . , M. Let 
us emphasize that we do not use the extension of fl in numerical algorithms and these 
extensions of f l {t) and fl are done in order to use the vector notation f(t) and f k without 
need to adjust length of these vectors as t and k grow. 

We assume that the A % in (13.151) are such that fl satisfy the following condition. 

Assumption 3.1 For ac£l there is a constant C > such that 

Eexp(c\fl\)<C (3.16) 

for alii = 0,...N, k = 0, . . . , M. 

This condition is satisfied by all sensible numerical schemes (i.e., sensible choices of 
A 1 in (I3.15P ) thanks to the uniform boundedness of <Ji(t,T,z) (see Assumption 2.1) and 
boundedness of the initial condition (see Assumption 2.3 and also the comment after 
(I3.13P ). In particular, it is satisfied by the weak Euler-type scheme (I5.3P we use in the 
algorithms in Section 
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We also require that the numerical method A3. 15)) for the SDEs fl3.10p - fl3.12p is of local 
weak order q + 1, i.e., that the following assumption holds. 

Assumption 3.2. We assume that the method (13.151) is such that for some positive 
constant C independent of A 

s s 

\E([[SP -l[5P)\ <Ch* + \ a = l,...,2g+l, (3.17) 

3=1 3=1 

2q+2 

E ]J \5p'\ < Ch q+1 , (3.18) 

3=1 

where 

<\i" ■= fUt + h)- x\ Sf := /;.,.(/ +h)- x\ 

and ft x (t + h) is the solution of the SDEs (13.101) with the initial condition x given at time 
t : fl x {t) = x\ and fl x (t + h) is the one-step approximation of (I3.10p found according to 
(jlriSD with fl x (t) = x\ 

Assumption 3.2 is similar to the one used in the standard theory of numerical inte- 
gration of SDEs in the weak sense (see, e.g. [IE])- As we will see in Section HI Assump- 
tions 2.1-2.3 and 3.1-3.2 guarantee weak convergence of the numerical method (I3.15P to 
the solution of the auxiliary system of SDEs (I3.10p with order h q . 

We note that C in (I3.17p - fl3.18p is independent of x while in the standard theory of 
numerical integration of SDEs one usually has C depending on x in such estimates (see 
[TBI P- 100]). In our case it is natural to put C independent of x since the coefficients 
of (13.101) and their derivatives are uniformly bounded (see Assumptions 2.1-2.2). We also 
emphasize that the constants C in (I3.17p - fl3.18l) are required not to depend on A. 

Remark 3.1 The numerical method <\3.15§ contains the approximation fl k of the forward 
rate f(tj.,Tt k ) (recall that > T^J which from the financial point of view does not exist 
unless tj. = Ti k . However, from both theoretical and numerical points of view, it is not 
prohibiting to consider the values f e k k which, as we will see later in Section l3~3\ is compu- 
tationally beneficial. We may interpret the points (tk,Ti k ) on our (t,T)-grid as fictitious 
nodes (see also Remark WJfy . 

The approximation (I3.15P of the infinite-dimension stochastic equation (I2.4p - fl2.5p has 
two discretization steps: T-step A (i.e., step in maturity time) and t-step h (i.e., step 
in time). We can say that the T-step A controls the error of approximating (I2.4p - fl2.5l) 
by (I3.10p - fl3.12l) while the t-step h controls the error of approximating (13. 100 - 03.120 by 
(I3.15p . We will later (see Remark 14. 3)) discuss how to choose A and h in practice. 

3.3 Approximation of the price of an interest rate contract 

In the previous section we introduced an approximation f l k of the solution to (I2.4p - fl2.5l) . 
Our aim is to evaluate the expectation (12.1 2ft - ( 12TT51) . i.e., 

F(t Q , f (■) ; s k , Si ) = Eexp{-Y{s k ))G {P{s k , Si )) . 
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To evaluate F, one has to compute Z(sk,si) from (I2.13P and Y(sk) from (I2.15p . In this 
section we construct numerical approximations for Z(sk, s») and Y(sk)- For clarity of the 
exposition, we assume in what follows that 

Sk = t* and Si = T* . 

We approximate the maturity time integral from ( 12. 1 5[) by a quadrature rule Sz(t*,T*, A) : 

Z(f,T*) = / f(t%u)duaS z (f,T*,A) = AY ( j j f(t*,T j ), (3.19) 

where the weights 7- are chosen so that the quadrature rule is of order p > 0, i.e., an 
inequality of the form ( 13. 6 P holds: 

(E [S z {t*,T*, A) - Z{t*,T*)} 2 ) 1/2 < CA P . (3.20) 

The assumption we made at the beginning of Section [3TT1 that there is a sufficient number 
of nodes 7] between t* and T* ensures that we can find a quadrature rule ( 13.191) satisfying 
(13 .2 Op . Some examples of quadratures Sz(t* ,T* , A) are given in Section |5j 

In general, T-discretization and t-discretization have different steps A and h, and 
approximate values of the short rate r(tk) = f(tk,tk) (cf. (12. 2ft ) are not directly available 
among fl which are defined on the (t, T)-grid. Then to numerically evaluate Y(t*), we 
need to construct an approximation of /(£&, £&) based on the values fl, i = £%>,... , N. To 
this end, let us first consider an approximation of the exact short rate r(t) = f(t,t) using 
the values of f(t,Ti), i = £(t), . . . , N. We recall that thanks to Assumptions 2.2-2.3 the 
solution f(t, T) of (I2.4[) -( |2T5T) is sufficiently smooth in the last argument. We approximate 
r(t) by 7r(t) as 

ir(t) = nit-Jit,^), i = £(t),...,£(t) + 9) (3.21) 

i{t*) 6 

e 

= ^2\i(t)f{t,T m+i ), te[t ,t% 

i=0 

where Aj(£) are coefficients independent of /, \\i(t) \ are bounded by a constant independent 
of A, 6 is a non- negative integer independent of t and A, and Xa 1S the indicator function 
of a set A. We choose the number 9 and the coefficients A,(t) so that the approximation 
( I3.2ip is of order p : 

[E [r(t) - n(t)} 2 ) 1/2 < CA P , p > 0. (3.22) 

The form of (I3.2ip covers both polynomial interpolation and extrapolation. For inter- 
polation, we approximate r(t) using the values f(t,Ti), i = £(t), . . . ,£(t) + 6. For ex- 
trapolation, the coefficient A (t) = and we approximate r(t) using the values /(£, Tj), 
i — @(t), . . . , g(t) + 9 — 1. Some particular examples of the approximation ir(t) are given 
in Section [5j Recall that in Section 13.11 we assumed that A is such that there is a suf- 
ficient number of nodes Tj between t* and T* which should, in particular, ensure that 
£{t) + 9<N. 
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Remark 3.2 We note that we need fictitious points (tk,Te k ) on our (t,T)-grid (see also 
Remark \3. 1\) for the interpolating form of (\3.21\i . The extrapolating form of ( \3.21\i does 
not need the fictitious points as it is sufficient to compute f k for i — g k , . . . , g k + 9 — 1, 
k = 0, . . . , M, all of which have the usual financial meaning. However, we reserve the 
possibility to use an interpolation (and, consequently, the fictitious points) for simulating 
short rates since interpolation is usually computationally preferable to extrapolation. 

Using the short rate approximation vr(s), we approximate the time integral in (12.131) 

as 

ft* t-t* _ ft* 

Y(t*) = / r(s)ds « / n(s)ds « Y{t*) := / n{s)ds, (3.23) 

J to J t() J to 

where tt(s) has the form of ir(s) from (13. 2 1 p but with f l (t) instead of f(t,Ti) : 

7f(s)=7r(s;/ i (s), i = e{s),...,e(s) + 6). 

We extend the system of SDEs (I3.10P by adding to it the auxiliary differential equation 

dY = tt(s; f(s), i = £(s), . . . , i(s) + G)ds, Y(t ) = 0. (3.24) 

Recall that n(s) for every s G [to, t*] is a linear combination of jf*(s), i = £(s), • • • , £{s) + 0. 

Let Y t>Xt y(s), s > t, be the solution of (I3.24[) with the initial condition Y t>X)y (t) = y and 
with p(s) = ft x (s) (recall that fl x (s) are defined in Assumption 3.2). We observe that 
under h < aA for some a > : 



t+ft " 

Y,Hs)ft:(s)Xs m ,T l+l) ds, (3.25) 

L=t[t) 1=0 

and for any positive integer m 

Y t , x , (t + h) <Ch m \l+ V \x l \ m I , (3-26) 




E 



where C > is a constant independent of A and x. We note that the condition h < aA 
guarantees that the number lit + h) — £(t) is independent of A, which ensures that the 
constant C in (13.261) is independent of A and the number of terms in the sum on the right- 
hand side of (13.261) is also independent of A. This will be essential for proving convergence 
Theorem 14.31 

Now we extend the fully discrete approximation (13.151) by adding to it an approxima- 
tion of (|g23|) : 

F = 0, Y k+l = Y k + A Y (t k -fl j=4,...,4+i + M), fc = 0,...,M, (3.27) 

where the form of A Y (i fc ; f k ,j = £ k , . . . , £ k+1 + 6;h) = A Y (t k ; h) depends on the form of 
7r(s) from (I3.2ip and the accuracy required. 

We replace Assumption 3.2 on the one-step approximation by the assumption which 
is applicable to the extended system (13.101) . (13.241) and the extended discretization (13.151) . 

dsizD. 
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Assumption 3.2'. Let h < a A for some a > 0. We assume that the method (13. 15)) . 
( I3.27P is such that for some positive constant C independent of A 



E ( 5Y m Yl Sf ij ~ SY m Yl Sf 

3=1 3=1 

m — 0, . . . , s, s — 1, . . . , 2q + 1 



£(t+h)+6 

< ( i + 

l=£(t) 



„l i m 



(3.28) 



E max 

0<m<2 ? +2,{ii,...i 29 +2-m}e{0,...,A r } 



2q+2— m 

5Y m [ [ Sf 

3=1 



-| 1/2 



< Ch q+1 1 



t{t+h)+e 

E 

i=£(t) 



J|2g+2 



(3.29) 



where 



sr = f tx { t + h)-x\ 6.r ./;.,(/ • /o 

= ^(t + /i) - y, 5Y = Y t>XtV (t + /i) - y, 

fl x (t + /i) and fl x {t + /i) are as in Assumption 3.2 and Y tjXjy (s), s >t, is the solution of 
(13. 241) with the initial condition Yt,x,y{t) — V, an d Y tyXyy {t + h) is its one-step approximation 
found according to (13.27)) with Y tjXjy (t) = y. 

Note that the constants C in ( 13 . 281) - (13 .291) do not depend on x, y, and A. The depen- 
dence of the estimates (13.28^ - (13. 29)) on x is consistent with (13.26)) . The condition h < a A 
in Assumption 3.2' is not restrictive from the practical point of view since we aim to be 
constructing efficient numerical algorithms for the HJM model by allowing bigger T-steps 
A without losing accuracy. We also note that this condition arises only when we need to 
approximate the short rate (see also Remark 14. 2D . 

Further, we make the following assumption. 

Assumption 3.3. For some c > and C > 



Eexp(c\Y k \) < C 



(3.30) 



for all k = 0, . . . , M. 

As a rule, the condition ( 13.301) immediately follows from Assumption 3.1 which is the 
case, e.g., for the algorithms presented in Section [51 

Based on (13. 19)) . (13.23)) and using (13 . 15)) . (13.271) . we arrive at the approximation F of 
F from (12~T2|) : 

F(t , f (•) ; t*,T*) w F{t , / ; t*,T*) = E exp(-Y M )G (P(t*,T*)) , (3.31) 
where Y M is from (13.271) ; 

P{f, T*) = exp (-S z (t*,T*, A)) , (3.32) 
Sz{t*,T*, A) is the quadrature rule of the form (13.191) with f(t*,Tj) replaced by f 3 M : 



N 



S z (t*,T*,A) = Aj2%f. 



■3 . 

Mi 



(3.33) 



3=Bm 
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and /o means the initial condition of (I3.10p . which is the N + 1- dimensional vector 
(fg,..-Jo N ) T = (fo(T ),...J (T N )) T . 

Finally, the expectation of the discounted payoff in (I3.3ip is approximated by the 
Monte Carlo method, i.e., 

F{t ,fo(-);t*,T*) « F(t ,fo;t*,T*) (3.34) 

L 



F(t J ;t*,T*) = i^exp(-y j S ) )G(P«(f,T*)) 



1=1 

where Y®, are computed using independent realizations fl , j =£&,..., iV, k = 
1 , . . . , M, of the random variables f 3 k . 

In f )3.34p the first approximate equality corresponds to the error of numerical inte- 
gration and the error in the second approximate equality comes from the Monte Carlo 
technique. The numerical integration error is analyzed in the next section. The Monte 
Carlo (i.e., statistical) error in (I3.34p is evaluated by 

_ _ [Var{eM-YM)G(P(t*,T*))}] l/2 

Pmc — c JJ/2 yo.oo) 



[Var {exp(-Y{t*))G (P(t*, T*))}} 



1/2 



L 1 / 2 

where, for example, the values c = 1,2,3 correspond to the fiducial probabilities 0.68, 
0.95, 0.997, respectively. The Monte Carlo error can be decreased by variance reduction 
techniques (see, e.g. [91 [TOj HU [20] and references therein). In this paper we deal with 
the numerical integration error and numerical algorithms which are effective with regard 
to (t, T)-discretization. 



4 Convergence theorems 



The aim of this section is to prove the convergence of the approximation F(to, fo]t*,T*) 
to F(t , f (•) ; t*, T*) as h and A ->■ 0. 

Denote by F(t , / ; t*,T*) the approximation of F(t , f (■) \ t*,T*) from fl2TT2|) result- 
ing from approximating the solution f(t,Ti) of f)2.4p -f l23|) by f l {t) from (I3.10p . i.e., 



F{t ,fo;t*,T*) = Eexp(-Y(t*))G (P(t*,T*) 



where 



P(t*,T*) = exp ( -S z {f,T*, A) 



(4.1) 



(4.2) 



Y(t) is from (ET23)1 and S z (t* } T*, A) is the quadrature rule of the form (13~T9|) with f(t*,Tj) 
replaced by f j (t*). 

The error R of weak approximation of F by F can be written as a sum of two con- 
tributing terms: 



R = F(t Jo(-);t*,T*)-F(t Jo;t*,T*) 
F(t ,f (-);t*,T*)-F(t Q J Q ;t*,T*) 

- R\ + i?2, 



(4.3) 



+ 



F(t ,f ;t*,T*)-F(t ,f ;t*,T*) 
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where R\ is the error due to T-discretization of (}2.4j) - (l2.5p and R 2 is the error due to 
t-discretization of (I3.10p - fl3.12p . The first error, Ri, is analyzed in Section [4.11 and the 
second error, R 2 , is analyzed in Section H~2l 

Note that in this section we shall use the letters K, C, and c to denote various constants 
which are independent of A and h. 



4.1 T-discretization error 



In this section we analyze the error of the finite-dimensional approximation f l3.10p - fl3.12p 
for the infinite-dimensional stochastic equation (12.4p - (l2.5p . The plan of this section is 
as follows. First, we prove (see Theorem I4.ip that the approximation f l3.10p -f l3.12p has 
mean-square convergence of order A p . This result plays an intermediate role for getting an 
estimate for the T-discretization error Ri but, at the same time, it has its own theoretical 
value. Based on Theorem 14. 1[ we prove (see Lemma [4. ip the mean-square convergence of 
Y from (I3.23P to Y from ( j2 . X3I) . Finally, in Theorem 14.21 we prove that the weak-sense 
error R 1 (see (JO}) of (jOTijl - fETIlZ)! is of order A p . 

Theorem 4.1 Suppose Assumptions 2.1-2.3 are satisfied. Then the approximation f l (t) 
from ( \3. 10§ - <\3. 12\ converges to f(t,Ti) from ( Jjg.^P - pTjp as A — > with the mean-square 
order p, i.e., 



E 



1/2 



< KA P , t e [t , t* A T ( 



(t+l)AJVj> 1 



0. 



(4.4) 



where K > is a constant independent of A, t, and i. 

Proof. Denote by p(t,Ti) the error of the approximation f l3.10p - p.12p : 



p^Ti) := f(t) - f(t,Ti), t e [h,t* AT (i+1)AJV ], z = 0,...,JV. 



Clearly, 



p(t Q ,Ti) =0. 

Due to Assumption 2.2, a(s, T, z) is globally Lipschitz in z whence 



laAs) - a(s,Ti 



and (cf. (J33D and ( B32D ) 

Si{s, Tj, A) — S/(s, Tj, A) 



a(s, Tt, r(s)) - a(s, T h f(s, 7-)) < K \p(s, T r 



K( S ,Ti) 

= A Yl ^k(s)(°k(s) -a(s,T k )) 

k=i(s) 

K{s,T t ) K( S ,Ti) 

< KA \a k (s)-a(s,T k )\<KA ^ \p(s,T k 

k=£(s) k=l(s) 



(4.5) 
(4.6) 

(4.7) 



(4.i 
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We have from (1XTU|) - (1XT2|) and (T2^j) -(l2~oD: 



p(t,Ti) 



to 



°l ( s )Si(s,T i} A) — a (s,Tj)/(s,Tj)J ds + / fc(s) - <r(s, T,)] 1 dW{s) 

ds 



to 



a, (s 



Si(s, Ti, A) — Sj(s, Tj, A) 
+ / a7(s) [S I (s,T i ,A)-I(s,T i )]ds 

+ [ [a l {s)-a{s,T i )} T I{s,T i )ds+ I [a t {s) - a(s,T,)] T dW(s). 

J to J t 

By Ito's formula, we obtain 



/ 2p( S ,T i )a7(s) b ) T 11 A)-S / ( S) T !) A) 

./to L 

+ / 2p(s,T l )a l (s) [S I (s,T i ,A)-I(s,T l )]ds 

J t 

+ / 2p(s, [afa) - a(s, Ti)] T I(s, T^ds 



ds 



to 



+ / [cri(s) - <t(s, T^)] t [ffj(s) - a(s, Ti)} ds 



to 



+ / 2p{s ) T l )[a l {s) - a{s,T t )}dW{s). 



to 



Then 



Ep\t,Ti) = 2 1 Ep(s,T,)aJ(s) Sj(s, T i} A) — Si(s, T i: A) ds (4.10) 



to 



+2 f Ep(s,Ti)drJ(s) [S I (s,T i ,A)-I(s,T i )]ds 

Jto 

+2 I Epis^ftiis) - a(s,Ti)] T I(s,Ti)ds 

Jt 

+ E [ai(s) - a(s, Ti)] T [ai(s) - a(s, T { )} ds. 



to 



(4.9) 



Using the boundedness of a(s,T,z) (see (12.61) ) and the inequality (14.80 . the first term on 
the right-hand side of (14.101) is estimated as 



2 f Ep(s, Ti)aJ( a ) Ti, A) - Sj(s, T u A)l ds 

Jto 1 J 



(4.11) 



< KA 



to 



K( S ,Ti) 

E\p(s,Ti)\ \p(*,T k )\ds 
k=e(s) 
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Using the boundedness of a(s, T, z), the inequality 2ab < a 2 + b 2 , and the condition (13. 6p 
for the quadrature rule 5*/, we obtain for the second term on right hand side of (14.101) : 

2/ Ep(s,T i )aJ(s)[S I (s,T i ,A)-I(s,T i )]ds (4.12) 

Jto 

< K I E\p(s,T i )\\S I (s,T i ,A)-I(s,T i )\ds 

Jto 

< K [ [Ep 2 (s,T l ) + E\S I (s,T l ,A)-I(s,T l )\ 2 ]ds 

Jto 

< K [ [Ep 2 {s,T l ) + A 2p ]ds. 

Jt () 

Using the inequality (14. Tf) and the boundedness of a(s, T, z), we get for the third term on 
the right-hand side of ( 14.101) : 

2/ Ep(s,Td[ai{s) - <r(s,Td] r I{s,fyds < K [ Ep 2 {s,T i )ds. (4.13) 

Jto Jt 

By the inequality ( 14 .7D . the fourth term on the right-hand of (14.101) is estimated as 

f E[d i {s)-a(s,T l )] T [d^-ait.T^dsKK I Ep 2 (s,T l )ds. (4.14) 

Jto Jto 

Substituting ( 14.1 1ft - ( I4T41) in (14.101) and using the inequality 2ab < a 2 + b 2 , we obtain 
Ep 2 ^,^) < K J\e P 2 (s } T 1 ) + AE 



n(s,Ti) 

\p(s,Ti)\ \p( s , T k 

k=i{s) 



A 2p } ds (4.15) 



< k[ { (l + A(K(s,T l )-e(s) + 2))Ep 2 (s,T t ) 

Jto 

K( S ,Ti) \ 

+A Ep 2 (s,T k ) + A 2p \ds 

k=£(s), k^i J 

,t [ ) 

< K lEp 2 (s,T l ) + A Yl Ep 2 (s,T k )\ds + KA 2p , 

[ k=e(s), k^i J 



t e [t ,t* AT { 



(i+l)AJVj, 



i — 0, .... N. 



We have used here that A(k(s,T;) - £(s) + 2) < T* - t Q . 

Introduce p M (t ) := m &x m <j< N Ep 2 (t,Tj), t <G [t Q ,t*]. Clearly (see (@2D)), p M (*o) = 0. 
Then we get from (14. 15ft : 



P M (t)<K f p M (s)ds + KA 2p , 

Jt 



whence ( 14.41) follows by the Gronwall inequality. Theorem 14.11 is proved. □ 
Using Theorem 14. 11 we prove the following lemma. 
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Lemma 4.1 Suppose Assumptions 2.1-2.3 are satisfied. The approximation Y(t) from 
(3.24) converges to Y(t) from (\2.13\i as A — > with the mean-square order p > 0, i.e., 



E 



Y{f)-Y(f) 



1/2 



< KA 1 



(4.16) 



where K > is a constant independent of A. 

Proof. Consider the error of the approximation (13. 24ft for (12.131) (see also (I3.23f) ) : 



Y{f) - Y(f 



t* i~t* 
f(s, s)ds — / n(s)ds. 

to J to 



(4.17) 



We rearrange the right-hand side of (14. 17ft to split this error into the error due to ap- 
proximation of the short rate r(t) = f(t, t) by 7t(t) and the error due to approximation of 



Y(t*)-Y(t*)= [ (f(s,s)-n(s))ds+ [ (tt(s) - tt (s)) ds. 

J t J t 



(4.18) 



Due to the condition (13.221) imposed on our choice of the approximation tt(s), we have 



(f(s,s)-n(s))ds S) j < KA 2p . 



(4.19) 



Recalling the form of the approximation tt(s) from (I3.2ip . we get 



E I I (tt(s) -n{s))ds 



E 



-1 2 



where A»(s) are bounded coefficients and the number 9 is independent of A. Then, using 
(j4~4j) . we obtain 

(4.20) 



£ (tt(s) - ?f(s)) ds^) < KA 2p . 



The relations (I4.18I) - (I4.20I) imply the required error estimate (14. 16ft . □ 

In the next theorem we obtain an estimate for the weak sense error Ri from (14. 3p . 

Theorem 4.2 Suppose Assumptions 2.1-2.3 are satisfied. Assume that the payoff func- 
tion G(z) satisfies the global Lipschitz condition (\2.11\) . Then the approximation F(to, fo, 
t*,T*) from (gjj converges to F(t J (-);t\T*) from p^ - pLXj with order 

p > 0, i.e., 

F(t J (-);t*,T*) - F(t J ;t*,T*)\< KA?, (4.21) 
where K > is a constant independent of A. 
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Proof. We have (see (EH2J), (ISTTHft - rtXTS]) and 

i?x = F(t , /o (•) ; ** , T*) - F(t , /„; t*,T*) 
= Eexp(-Y(t*))G(P(t*,T*))-Eexp(-Y(t*))G (P(t*,T* 

= E exp (-y(t*)) - exp (~y(t*)) G(P(t*,T*)) (4.22) 

+£ [G (P(t*,T*)) -G(P(t*,T*))] exp(-Y(t*)). 



Consider the first term on the right-hand side of (14.220 . By the mean value theorem, 
we get 

exp (-Y(t*)) - exp (-Y(t*)) = (Y(t*) - Y(t*)) exp(tf), (4.23) 



where $ is a point between —Y(t*) and —Y(t*). Due to the global Lipschitz condition 
(12. lip imposed on G(z), we have (recall that Sz(t*,T*, A) is the quadrature rule of the 
form (13411 with /(f ,T;) replaced by /*(**) : 



|G(P(t*,T*))| < KP(t*,T*) = K exp (—S z (t*,T*,A) 



(4.24) 



A' 



Kexp -A ^ 7, 



Using (14.231) . (14.241) . and the Cauchy-Bunyakovsky inequality twice, we obtain 



E 



< K 



exp(-y(t*)) -exp(-y(0) G(P(t*,T*)) 

1/2 



E(y(f) - F(r)) J [£exp(4tf)] 



,1/4 



iV 



£exp -4 A J2 7i W 



(4.25) 

1/4 



Thanks to (12. 8p and (I3.13p . exponential moments of —Y(t*) and —Y(t*) are bounded and, 
consequently, for some K > we get £"exp(4^) < K. Due to (I3.13p . we also have 



N 



Eexp -4A 7i 



)-« 


exp | 







iV 



exp -4T ]T 7, /'(**) 



3=0 M 



1/N 



< 



1 ^ 

-^£exp(-4T7^(t*))<K 



o=e M 



Then (I4.25P together with (14.160 implies 



\E 



exp (-y(r)) - exp(-y(f )) G(P(t*,T*))\ < KA P . 



(4.26) 



Let us now consider the second term on the right-hand side of (14.220 . Due to the 
global Lipschitz condition (12.110 imposed on G(z), we have 



G(P(t*,T*))-G(P(t*,T*)) 



<K- 



P(t*,T*)-P(t*,T*) 



(4.27) 
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Further, by the mean value theorem, we get 

P(t*,T*)-P(t*,T*) = exp (-Z (f* , T*)) — exp (—S z (f, T*, A) 

/ S z (t*,T*,A)-Z(t*,T*))exp {<&) , 



(4.28) 



where is between -S z (t*, T*, A) and -Z (t*,T*). Using (OTjl . (fOSjl . and the Cauchy- 
Bunyakovsky inequality twice, we obtain 



G (P(r, T*)) - G(P(t*,T*)) exp(-Y(t*)) 



(4.29) 



< 



£ (S z (t*,T*, A) -Z(t* ,T*; 



1/2 



[Eexp(-4F(r))] 1/4 [£exp(4$)] 



,1/4 



It is clear that (12.81) and (13.131) imply boundedness of the exponential moments present 
in the right-hand side of (14.291) and, hence, 



E 

< K 



G (P(t*,T*)) - G \P(t* , T*)JJ exp(-y (t* )) 

, . 9.1 1/2 

£ (5 z (t*,T*, A) -Z(t*,T*] 



(4.30) 



We have 



£ (S z (t*,T*, A) -Z(t*,T*) 



(4.31) 



£ (5 z (t* , T*, A) - S z (t*,T*, A) + S z (t*,T*, A) - Z (t*, T*; 



< 2£ 



^(r,T*,A)-^(r,T*,A) 



+ 2£ 



S z (t*,T*, A) -Z(t*,T*) 



n 2 



Due to the condition (13.201) imposed on the quadrature rule S z (t*,T* } A), the second term 
on the right-hand side of (14.311) is bounded from above by KA 2p . Using (I4.4p . we obtain 
for the first term on the right-hand side of (14.311) (cf. (13.191) ): 



IE 



S z (t*,T*,A)-S z (t*,T*,A) 



2A 2 E 



N 



N 

< KA J2E[f(t*)-f(t*,T 3 



3=8 At 



< KA(N - g M + l)A 2p < KA 2p . 



Hence 



E (§ z (t*,T*,A) - Z(t*,T*)) <KA 2p . 



(4.32) 



The required estimate pill follows from (14^22]) . fT4T26]) . (00]) . and (1432]) . Theo- 
rem 14.21 is proved. □ 
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4.2 ^-discretization error 



In this section we analyze the error R2 (see (I4.3P ) due to t-discretization of (I3.10p -f l3.12p : 

R 2 = F(t ,f ;t*,T*)-F(t ,f ;t*,T*). 

Then combining its estimate with the estimate (I4.2ip for R\ from Theorem 14.21 we prove 
convergence of the weak approximation F to F (see (14. 3p ). In the analysis of R2 the key 
is to show that convergence of F(t , f ;t* ,T*) to F(t , f ;t* ,T*) is uniform in A, which 
is the reason why we cannot just apply here the standard results of weak convergence of 
numerical methods for SDEs (see, e.g. [IB])- The convergence theorem is proved under 
the assumption that the pay-off function G(z) in (I2.12p is sufficiently smooth. At the end 
of this section we also discuss how this assumption can be relaxed. 

To prove the convergence theorem (Theorem 14.31) of F(t , f ; t*,T*) to F(t , f ; t*,T*), 
we need the following technical lemma. We will use the multi-index notation: 

i = (z , . . . ,In) 

with ij being nonnegative integers, |i| = i + ■ ■ ■ + i^, and i! = i \ ■ ■ ■ ij^l. 
Lemma 4.2 Let A m be the m th -order operator 

gm 



A m = A? = V u 1 

^ (f)<r-0 
\i\=m 



^ {dx°Y° ■ ■ ■ {dx N f N 



with some fi 1 . Suppose Assumptions 2.1 and 2.2 are satisfied. Assume that the payoff 
function G(z) has m* bounded derivatives. Then for m > up to the order m* 

A m F(t,x;t*,T*)| < ^ Max exp(cA|x|), (4.33) 

where K > and c> do not depend on A and x G M^" 1 " 1 , and fJ^Max := max |i|=m l/^'l- 

Remark 4.1 To help with intuitive understanding of this lemma, we remark that A m F 
can be viewed as a Frechet derivative of the option price with respect to the discretized 
initial forward rate curve. 

Proof of Lemma 14.21 Recall the notation: ff x (s), s > t, is the solution of the system 
of SDEs (I3.10p - fl3.12p with the initial condition at t > t : fl x {t) = xK We introduce a 
more detailed notation for Sz{t*,T*, A) (cf. (I3.19P ): 

N 

S z (t,x;t*,T*,A)=Aj2 JjfUt*), (4-34) 

3=Qm 

which we can present as (cf. f)3.10p ) 

N 

'<■./ 



S z (t,x;t*,T*,A) = A £ V 

0=Qm 

N r r f ,/ 



+A E^' / <T T (s,T j Jl x (s))S I (t,x;s,T j ,A)ds + £ a 1 (s,T„ fl x (s))dW(s) 



3= Qm 
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where (cf. <KT2\> ) 



K(s,Tj) 

S I (t,x;s,T j ,A) = A £ 7i(*M«, T h 

i=e(s) 

Then, thanks to Assumption 2.1, we obtain for any positive integer m 

EP m {t*,T*) = Eexp (-mS z {t,x;t*,T*,A) 



(4.35) 



(4.36) 



N 



= Eexp J —mA 7 J x- ? 
\ J=s M 

" mA Z)Ti[r cx T (s,T j Jl x (s))S I (t, x- s,T j ,A)ds+ f a T (s,T v f~l x (s))dW(s) 

/ ,v 
< i^exp J cA | re 



JV 



3=Qm 



1=Qm 



where K > and c > do not depend on A. 

Further, recall that Y ttXty (s), s > t, is the solution of (I3.24p with the initial condition 
Yt,x, y (t) = y and with f l (s) = /^(s), i.e., 



Y t , x ,y(s) 



= y + J ft(s)ds 

y+ r^s'-jUs'), i = t(s'), 



(4.37) 



(s') + 9)ds' 



l=t(t) i=Q 



5) z+e 



^EE 

|=£(t) m=( 



sAT i+1 



iVT, 



sAT !+1 



iVT, 



^ m -i(s')fZ(s')ds', t>t Q , s>t, 



where (cf. (13.211) ) 6 and A m _;(s') depend on our choice of the accuracy order of short rate 
approximation, and 6 does not depend on A, and |A m _{(s')| are bounded by a constant 
independent of A. 

We also see that Y tjXjV (s) = y + Y ttXt0 (s). Using (I4.37p . (I3.10p . and Assumption 2.1, one 
can show that for any x > 



E 



ex P (x|y^ (r)|) 



.Eexp I x 



l=e(t) m=l JtyTl 



X m -i(s')fZ(s')ds 



(4.38) 



i(t*)+e 



< Kexp cA ^ 

i=i(t) 



\x 



where K > and c > do not depend on A. 
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Using smoothness of G(z), we obtain 



A m F(t,x;t*,T*) = EA m ex P (-Y tiXfi (t*))G(P(t*,T*)) 
E V // : -exp(-Y tx0 (t*))G(P(t*,T*)) 



(4.39) 



\i\=m 



(dx y° ■ ■ ■ (dx N ) 



m m—k t m 



x— G (P(f,T*)J P p (t*,T*) 



N 



fc. 



X 



n 

n=l 



—S z (t,x;f,T*,A), 



where C(a,/3,ji, . . . , j ktJ h, ■ ■ ■ ,l n *) are constants independent of N; j k = ^2 r=1 j r , l n = 
Ylr=i hi the sum Yli~j k +T n = m * s taken over all positive integers ji, . . . , and l±, . . . , l nt 
such that j k < j k+1 , k — 1, . . . , fc* — 1, Z n < / n+ i, n = 1, . . . , n* - 1, and j kt + [„„ = to; 
and in the right-hand side the multi-index i at /x 1 corresponds to the values taken by 
k, ■ ■ ■ , i] kt , r l7 . . . , . 
We have (cf. (Oil ): 



5' 



A? 



5' 



7c 



Q=Bm 



dx {l ■ ■ ■ dx il 



/&(0- (4-40) 



Using the Cauchy-Bunyakovsky inequality, the assumed boundedness of derivatives of 
G(z), and the inequalities fl4T3oD and (jQ5]l . we obtain from flCTD - flQOj) : 



2V 



|A w F(t,x;f ,T*)| < ifexp cA ^ 



(4.41) 



m m—k t 



x EE E I* 



A 



e /'n^7 



il}'"^Ti ^1,— ,rr ^ =0 fc=l 



A 



X 



n 

n=l V g=e M 



7 



5' 



9 <9x ■ ■ ■ dx 



-fun 



dx 1 ** 

2\ 1/2 



where K > and c are independent of A and x. Then, to complete the proof of this 
lemma, it is sufficient to show that for any < fc* < m and < n* < m — fc*, any 
combinations of ji, . . . ,j k * and li, . . . , l n * satisfying j ksf + l n * = m, and any combination 
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of q x , . . . , q n * with Q M < qi < N : 

N 

e «'n 



E 



k, 



Q3u 



Yt,x,o(t*) 



(4.42) 



n 

n=l 



5'' 



where K > is independent of A and x. 

We can obtain the following SDEs (see (13.241) ): 



d 



dx 11 ■ ■ ■ dx 



dx n ■ ■ ■ dx 
and (see (ETTDjl ) 

d 



j-Yt t x,o{t) 



t(t*) i+e 

l=£(t) r=l 



0. 



dx ix ■ ■ ■ dx 1 



-fL(s)ds, 



d 1 



l l l—l l—n* 



dx ri ■ ■ ■ dx n 



a=0 /3=0 n,=0r,=l +p T<t =i 



K(s,T q ) 



X 



r/ a ~ d@ 

E A E 7. ^ T (^^^(«))^ T (^^ ) /^(«)) 

{fc 1 ,...,fc i }={ri,...,r ! } t>=£(s) 



n=1 <9x ■ ■ ■ dx 



i i 



_~[ 9a; i+'n. +Pr-i . . ■ dx ln 
d a 



a=l n.=l l=z 



n 



{fei,...,fc;}={n,...,r i } 



:*>(8,T q ,fL( 8 )) 



n=i dx^ 1 "- 1 ■ ■ ■ dx k 

Ql 



-fL(s) dW(s) 



dx ri ■ ■ ■ dx n 



flx( s ) — Xl=ll 



where C(a, (3, n*, r*) and C(a, n*) are constants independent of N, and ^2^ kl ki}={n n} 
means summation over all possible recombinations {ki, . . . , ki} of r%, . . . , 77 (note that the 
number of terms in this sum depends on I but not on N) . 

To obtain (14.421) . we first consider the case m — 1 for which it is sufficient to get 



an estimate for E 



J2iLo ft,x( s ) ■ To this end, introduce the process Ct,x( s ) 



(C u (s),...,C JV (s)) T with C(s) := EiU^mUs), s > t, which satisfies the following 
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system of SDEs 

d( j = ^(s^Jl^s)) ■ Si(t, x; s, Tj, A) • C j ds 

k(s,T) 

+° T {s ) T v ~flM) • A £ ^{s^Jlis)) .( l ds + ^is^flMX 3 dW{s), 

l=l(s) 

( j (t) /<'• ./ o v. 

Then using Ito's formula and Assumptions 2.1 and 2.2, we obtain after some straightfor- 
ward calculations: 



rs PS 

E[(IM 2 <K[^] 2 + K E[C(s')] 2 ds' + K Aj2 E [( l ( s ')] 2 ds'. 

Jt Jt '-t(s') 



Let £{s) := maxo<j<jv E \Ci x {s)\ ■ Then 

£{s)<K^ Max + K / £(,')<!,' 



where K > does not depend on A and x. Hence, by Gronwall's inequality 

£(s)<Kii 2 Max , t<s<t*. (4.43) 
Next, we consider the (N + l) 2 -dimensional process 



U,i2=0 



2 

Using the same recipe as in the case of estimating max <j<Ar E [Ct,x( s )] > we § e t that 

max £ [Ci 1 ^ 2 (a)] ' < ^L, (4-44) 

0<j<N 

where K > does not depend on A and x. Using (j4.44j) and repeating the same recipe 
again in the case of the processes 



02 



2 



dx %1 dx 12 ' 

»1,12=0 

- 

11,12=0 



we obtain 



max E [ 2 C{ x (s)(s)] 2 < K^ Max , (4.45) 
max E[rf tjX (s)] 2 < K^ Max , (4.46) 
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where K > does not depend on A and x. Using (14.451) . it is not difficult to get that for 
the process 



02 



ii,i2=0 

the following estimate also holds 

E[ 2%x {s)] 2 <K^ Max . (4.47) 

It is clear that ( 1Q4l) -f Q71) are sufficient for proving (02"]) with m = 2. To show ( 14421) 
for m = 3, we need to obtain estimates for the second moments of the processes 

U,i2,«3=0 

iV 02 o 

«1,*2,*3=0 



3 



= E f'™ ^dx*dx* & {s) > j= °'---> N ' 

ili*2,*3=0 



and 



«1,*2,«3=0 

iV f) 2 f) 

*Ais) = E ^ 1,l2 "^^^ (s) ^- o(s) ' J = '---' iV ' 

«1,*2,«3=0 

ii,«2,i3=0 

- «9 3 

11,12,13=0 

which can be done using the same recipe but with more laborious calculations. In the case 
of an arbitrary m one need to consider processes ( n, '"' Jm (s), zC'tx '" Jm ~ :L (s), • • ■ , m(t,x( s )i 
V^x '"'* ,m ~ 1 (s), m _i j i^ x (s), . . . , mVtx( s ) defined in the same fashion as we did in the cases 
m = 2 and 3. It is not difficult to see that employing the same recipe maxima of their 
second moments will be again bounded by Kjj? Max1 from which (I4.42p follows for an 
arbitrary m. 

The required inequality f 14 . 3 3 1) follows from (14.41 j) and (14 .42 p . Lemma 14.21 is proved. 

□ 

Using Lemma S21 we now prove convergence of F(to, fat* , T*) to F(t , fat* , T*) in 
the case of smooth payoffs G. 

Theorem 4.3 Let h < a A for some a > 0. Suppose Assumptions 2.1-2.3 and Assump- 
tions 3.1, 3.2', and 3.3 are satisfied. Assume that the payoff function G(z) has bounded 
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derivatives up to a sufficiently high order. Then the approximation F(to, /q; t*,T*) defined 
by (fOIjl - pjgj) . fTQ3j) . ( g^ ) conwerc/es to F(t J ;t\T*) from with order q > 0, 

z.e., 



F(t Q J ;t*,T*)-F(t J ;t*,T* 
where K > is a constant independent of h and A. 



< Kh q , 



(4.48) 



Proof. Using the standard technique (see [HI p. 100]), we can write the difference R2 
in the form 



R 2 = F(t ,f ;t*,T*)-F{t J ;t*,T*) 

= E exp(-Yt ,fo,o(tM))G(eM-Sz(t , /o; t*,T*, A))) 
-Eexp(-Y M )G(eM-Sz(t*,T*,A))) 



(4.49) 



M-l 



£ £ [exp(-^ )/l) y i (* M ))G?(exp(-5z(* i ,/ i ;t*,r, A))) 

i=0 

-exp(-y ti+1)/f+1 ^ +1 (t M ))G(exp(-^(t m ,/ i+1 ;t*,T*,A))) 



M-l 



£E{exp(-n iJ ^ i (t, + 0)^[ex P (-^ +i4j _ (ti+i)i0 (t M )) 



i=0 



xG(exp(-^(t, +1 ,4 )/i (t i+1 );t*,T*,A)))|/ ti)/i (t l . +1 )] 

-exp(-F m )E[ex P (-l> i+lJi+li0 (tM))G(exp(-^(t m ,/ i+ i;r,T*,A)))|/ m ]} 



M-l 



££{exp(-:P tiJ ^(^ 

i=0 
M-l 

£ £exp(-F,)£ [exp(-y t4i/ii0 (t i+1 ))F(t i+1 ,/ ti , /l (^ +1 );r,T*) 



i=0 



-exp(-y ti ,. )0 (t i+1 )) J P(t i+1 ,4 i/i (t m );t*,T*)|/ i 



M-l 



i=0 



where 



p(t,x) = E[exp(-Y t<xfi (t + h))F(t + h,fl x (t + h) ] t*,T*) 
- exp(-Y t:Xfi (t + h))F{t + h, f t , x (t + h); f,T*)]. 



(4.50) 



Now we write the Taylor expansion of the terms under expectation in (I4.50p in powers of 
5Y = —Y t:XjQ (t + h) and 8f % = f\ x {t + h) — x l and in powers of SY = —Y tiXtQ (t + h) and 
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Sftx = ftx(t + h) — x\ As a result, we obtain 



exp(-y t , Xi0 (t + h))F(t + h, f t , x (t + h);t*,T* 
F(t + h,x;t*,T*) 

29+1 ^ d \i\ 



(4.51) 



+ ^ i!Jfc! 



|i|+fc=l 



+ iijfcj 



(dx°) i0 ■ ■ ■ (dx N ) iN 



i0 / ~j\A %N ~ i 
- ~N\ 5 yk 



\i\+k=2q+2 



(dx°y° ■ ■ ■ (dx N Y N 
x exp{-m^, {t + h)) x (5f°y° ■ ■ ■ (5f N J N 6Y k 



F(t + h,x;t*,T*)(6f) •••(5/ 
F{t + h,x + x{Mt + h)-x);t*,T 



where x an d are from [0,1]. 
Further, 



eM-Y t , x ,o{t + h))F{t + h, f t , x {t + h); t*,T*) 
F(t + h,x;t*,T*) 



(4.52) 



2q+l 



+ N 5 =1 m (dxT ■ ■ ■ (dx N ) 



+ ... i!fc! f,9x 



|i|+fc=2(j+2 



(da; )* • ■ ■ (dx N ) lN 



F(t + h, x- £*, T*) (5f°) io ■ ■ ■ (5f N ) iN 5Y k 
F{t + h,x + xiftAt + h)~ x);t*,T*) 



x exp(-9Y t>Xi0 (t + h)) {SfT ■ ■ ■ {5f N Y N 8Y k , 

with x an d $ being from [0, 1]. 

It is not difficult to check (see also (I3.26P ) that under the assumed condition h < a A, 
a > 0, the following inequality holds: 



E max 

0<m<2g+2,{ii,...i 29+ 2-m}e{0,... i 7V} 



2q+2-m 

5Y m Yl 5 P 



1/2 



l(t+h)+6 



< 



Ch q+1 1 + 



\x 



l\2q+2 



l=l{t) 



(4.53) 



where C > is independent of A. We note that the number of components x l appearing 
in the right-hand side of (14.531) is not larger than 1 + 8, which does not depend on A. 

Using Lemma I4T2| the inequalities (I4.53p . (I4.38P and ( 13 . 1 3j) . and the Cauchy-Bunyakovsky 
inequality, we obtain 



9 \i\ 



|i|+fc=2g+2 



\k\ (dx°y° ■ ■ ■ (dx N f N 



-F(t + h,x + xifUt + x);t*,T*) (4.54) 



x exp(-6Y t , xfi (t + h)) (5~f°y ■ ■ ■ (sf N J N 5Y k \ 
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2q+2 

< E exp(\Y ttXfi (t + h)\) x J2 

am 



fc=0 



|i|=2<j+2-fc 



'o 



Sf N )" 5Y" 



IN 



(dx ) 10 ■■■(dx 



F(t + h,x + xiftA* + h)-x)- f,T*] 



< KE 



expfli^oCt + JOl) m ax \(5f° 

|i|+fc=2g+2 V 



5f N 5Y k 



x exp(cA|x + xUt )X (t + h)~ x)\) 

£(t+h)+6 

< Kexp{cA\x\) ( 1+ l x '! 29+2 

l=£{t) 



where K > and c > independent of A, h, and x. 

Analogously, using Lemma I4.2[ the inequality (I3.29|) from Assumption 3.2', and As- 
sumptions 3.1 and 3.3, we get 



9 \i\ 



F(t + h,x + x{ft,S + h)-x); t*,T*) (4.55) 



x eM-8Y t>x ,a(t + h)) {5f°y° ■ ■ ■ {Sf N y 

(t(t+h)+e \ 
1+ \A 2q+2 \h q+1 . 
l=£(t) J 



We obtain from fl430|) - fr4T52|) and (jO^I . ffl~55D : 



2q 



\p(t,x)\ < 



k=0 



2q+l-k 



5 m 



^ ^ k (dx ) 10 ■ ■ ■ (dx N ) 



IN 



F(t + h,x;t*,T*) 



(4.56) 



e(t+h)+e 

+Kexp(cA|x|) ( 1+ \ xl \ 29+2 1 



with > and c > independent of A and 
1 



E (spy ■ ■ ■ (df N y N 5Y k - E (5f°Y° ■ ■ ■ (5f N ) iN SY l 



Applying Lemma 14.21 and using the inequality f!3.28j) from Assumption 3.2', we obtain 
from fl4~56l) : 

e{t+h)+e 



\p{t,x)\ < ifexp(cA|x|) j 1 + \ xl \ 2q+2 ) h9+1 > 

i=l{t) 



(4.57) 



where K > and c > do not depend on A and x. 

Substituting (14.571) in (I4.49[) and using Assumptions 3.1 and 3.3 and the Cauchy- 
Bunyakovsky inequality, we arrive at the required (I4.48p . Theorem 14.31 is proved. □ 
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Remark 4.2 As it follows from the proof, in Theorem \4-3\ the condition h < a A is used 
only for estimating the parts of the error involving the approximate discounting factor 
Yt,x,o( s )- If pricing an interest rate derivative does not require a discounting factor (e.g., 
when one uses the forward measure pricing, cf. Remark 1 2. 2\) then a theorem analogous 
to Theorem \4-3\ can be proved under Assumptions 2.1- 2.3 and Assumptions 3.1 and 3.2 
without the restriction on h. 



Theorems 14.21 and 14.31 imply the following result. 



Theorem 4.4 Under the conditions of Theorems \4 ■ 1\ and \4 ■ 3[ the approximation F(to, /o; 



t*,T*) defined bu ffXI7j) - fTQgj) . (fX75j) . pi converges to F(t Q , f (•) ; t*, T*) from pOD 



<\2.15T\ . t\2.4\ with order p > in A and with order q > in h, i.e., 

\F(t , f (•) ; t*, T*) - F(t , /„; t*, T*) \ < K(A P + h% (4.58) 
where K > is a constant independent of A and h. 

Remark 4.3 (Relationship between A and h) A higher order p, i.e., a higher order of an 
approximation F(to, fo; t*, T*) of F(to, fo (■) ;t* ,T*), can be achieved by using a higher- 
order quadrature rules in (\3.9i) and (\3.19i) and higher- order interpolation or extrapolation 
in {\3.21\j . For this purpose, we can use a large arsenal of effective quadrature rules and 
interpolation/extrapolations methods from the deterministic numerical analysis which are 
directly applicable here (see Section^. To achieve a higher order q, we need a higher-order 
weak-sense numerical scheme for (\3. lLk - <\3. 12\\ . As it is known (see, e.g. fJEf), this is a 
harder task, and, due to complexity of stochastic schemes, one usually restricts themselves 
to using weak methods of orders 1 or 2. As a result, in practice we will takep > q. Then, to 
balance the two errors in <\4-58 ), we choose A = ah q l p for some a > to obtain the overall 



error to be of order 0{h q ). In other words, by increasing the order p we can take larger 
T -discretization steps A and, consequently, significantly improve computational efficiency 
of HJM simulation which, in particular, is illustrated in our numerical experiments in 
Section [7| 

According to the motivation examples considered in Section [21 the payoff G(z) is 
usually globally Lipschitz (see (12. f f h ) but not sufficiently smooth function as it is required 
in Theorem 14.31 and, consequently, in Theorem 14.41 Let us discuss two ways how one can 
deal with this theoretical difficulty. 

First, as it was noted in, e.g. |19j . we can approximate the payoff function G(z) by a 
smooth function G(z). Denote by e an error of this approximation. The proposed numeri- 
cal method can be applied to the smooth approximating function G(z) and Theorems 14.31 
and 14.41 remain valid for F with G instead of G. In this case, the overall error in evaluating 
the price of an interest rate contract consists of the numerical integration errors estimated 
in Theorem 14.41 and the error e of the smoothening of G. 

Second, one can exploit the result of [Tj which in application to our problem means 
that if the transition Markov function for the process f(t) is sufficiently smooth and fl 
is simulated by the strong Euler scheme then F(t , f ;t*,T*) converges to F(t , f ;t*,T*) 
with order one in h even for nonsmooth G. 

We remark that the computational practice (see our numerical experiments in Sec- 
tion [7]) suggests that the error estimates of Theorems 14.31 and 14.41 are valid for the weak 
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Euler-type scheme (see (15.31) below) in the case of nonsmooth G(z). Further, it is natural 
to expect that for higher-order weak schemes the error estimates of Theorems 14.31 and 14.41 
are also valid for nonsmooth payoffs G(z). We note that to answer on these theoretical 
questions related to nonsmoothness of G(z) further development of the general theory 
of numerical integration of ordinary SDEs is required which is outside the scope of the 
present paper. 



5 Numerical algorithms 

In this section we provide some particular examples of the generic numerical method 
introduced in Section [3J For simplicity of the presentation, we restrict ourselves in this 
section to the case of T-step being not larger than the t-step, i.e., 

h<A, (5.1) 

which is a stronger condition than the one assumed in Theorems I4.3l and l4.4l : h < a A, a > 
0. This requirement is not particular restricting since our aim is to construct efficient 
algorithms for the HJM model by allowing bigger T-steps A without losing accuracy as it 
was discussed in the Introduction and Remark 14.31 We note that there is no difficulty in 
constructing algorithms imposing h < a A for some a > instead of (15.11) . The condition 
(15. ip ensures that there cannot be more than one node Tj in any interval [t k ,t k+ i) hence 
there are only two cases possible: either 4+i = 4 or 4+i = 4 + 1. This is used in 
constructing numerical algorithms of this section. 

We need the following new notation in this section: 

A itk :=T t -t k (5.2) 

and 

4 + t k+ i 
tfc+i/2 — ^ ' 

In the paper we limit the illustration (see also Remark 15 .3 j) of the generic numerical 
method from Section [3] to considering only the weak Euler-type scheme (i.e., with q — 1) 
as a numerical approximation of the SDEs (I3.10p . (I3.24p . i.e., as an approximation of the 
t-dynamics. In this case the extended discretization (I3.15p . (I3.27P takes the form 

/* = /», i = 0,...,N, F = 0, (5.3) 
fUi = f l k + E <M*k)S/,(*fc, Tn A, h) + h 1 ' 2 £ rr, ; (/,K ; .,.,- 

3=1 3=1 

i = 4+i, ■ • -,N, 

Y k+1 = Y k + A Y (t k ; ]{, j = 4, • • • , £(t*) + 8; h), k = 0, . . . , M - 1, 
where £ ■ k+l are independent random variables distributed by the law P(£ = ±1) = 1/2, 
(&i,i(tk), ■ ■ ■ , o"i,d(4)) T = (cri(tk,T h fl), . . . ,a d {t k ,T h f k )) J , 

(tk, Ti] A, h) depends on our choice of the quadrature rule (I3.9p . and A Y is as in (13.271) 
and depends on the choice of approximation for the short rate (13.211) . 

In the remaining part of this section, we give three algorithms based on rectangle (p = 
1), trapezoid (p = 2), and Simpson (p = 4) quadrature rules SrAtk, T, A) accompanied by 
short rate approximations of the corresponding orders. In all these cases it is not difficult 
to check that (15.31) satisfies Assumption 3.1 and that Y k satisfy Assumption 3.3. 
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5.1 Algorithm of order 0(A + h) 

The application of the composite rectangle rule to approximate the integrals Ij(tk,Ti) in 
( 1511) and Z(t*, T*) in flUHJ) yields 

Sj.(4,T 4+1 ;A,/i) = /iA i k +i,kCe k+1 ,j{tk), (5-4) 
IiK k ' " k+v ' J I S J .(t fc ,r w ;A,/i)+AAa ffjb+1>j (t fc ) > otherwise, 

i 

Si.(t k ,Ti;A,h) = §i.(t k ,T ek+1 ;A,h) + hA ^2 °™.j{tk), « = + 1, • • • , 

m =e fc+ i+i 

^(r,T*,A) = /l-A, MiM + A £ (5.5) 

By straightforward calculations one can show that the used rectangle rule satisfies the 
order conditions (I3.6P and (I3.20p with p — 1. We pay attention that we incorporated two 
cases in (j5.4p : when £ k+1 = £ k and hence T lk+1 < t k and when (see also fl5.ll) ) £ k+ \ = + 1 
and hence Ti k+1 > t k . 

We use the piecewise approximation of the short rate (cf. (13.2 lj) ) : 

t{t*) 

n(t) = J2f(t,T l ) Xte[TitTi+i) , te[t ,n (5.6) 

1=0 

The approximation (15.61) obviously satisfies the order condition (I3.22p with p = 1. To 
satisfy Assumption 3.2' with q = 1, we, in particular, need to approximate the integral 
Y t , x , (t + h) in (13725]) by Y t)X)0 (t + h) from (EOT)) with local order 0(/i 2 ). In the case of flgjlfl 

the coefficient in the right-hand side of (13.241) ir(s;x l , i = £(s)) = / J^XsefT; ,t m ) = x ^ 

1=0 

is only piece- wise smooth. Further, according to the condition (15. ip . we can have two 
cases: either an open interval (t k ,t k+ i) does not contain any node Tj of the T-grid or it 
contains a single node T g . In the former case we can approximate the integral Y tiXi0 (t + h) 
in (13.251) by the left rectangle rule and we have A Y (t k ;fi, j = £ k ,£ k+ i;h) = hf h h with 
the local error of order 0(h 2 ) as needed. In the second case to achieve the local error 
of order 0{h 2 ) despite lack of smoothness of 7r(s;x l , i = £(s)), we split the integral 
Yt,x,o(t + h) = Y ttXfi (T £k+1 ) + Yr gk j(T ih+1 ),o(t + h) and approximate the first integral by the 
left-rectangle rule and the second by the right-rectangle rule: A Y (t k ; f£, j — £ k , £ k +i, h) = 
Ae k+1 ,kft - A^+i^+i/^ 1 - Tnus > 

A Y {t k ; f k ,j = 4,4+H h) = (hA A w ) ft - (0 A A 4+ljfc+1 ) (5.7) 

We note that despite the use of f k ^ in the right-hand side of (15.71) the method does not 
require to resolve any implicitness. 

Assumption 3.2' with q = 1 can be checked for the scheme (15. 3p . (15. 4p . (15.71) following 
the standard, routine way (see, e.g. [TSJ Chap. 2]). 
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The algorithm based on ( 15. 3 p and (15. 4p . ( 15. 5p . (15. 7p . we will call Algorithm 5.1 for the 
option price ( I2.12p -( l2.15p . According to Theorem 14.41 this algorithm is of order 0(A + h), 
which under the condition (see (15. ip ) A = ah, a > 1, resulting in 0(h). We also note that 
in the case A = h the short rate is readily available on the grid and its approximation 
is not needed. Algorithm 5.1 with A = h is analogous to the numerical methods for 
the H JM model considered in [T2~| fT5| [9] . As it is shown in our numerical experiments 
(see Section [7]), Algorithm 5.1 is less efficient than the new algorithms (Algorithms 5.2 
and 5.3) which we propose in the next two sections. 



Remark 5.1 If we replace A Y in (|5. ?p by 



A Y (t k ;h) = hf k * (5.8) 

then Assumption 3.2' with q = 1 is not satisfied and we cannot guarantee closeness of 
Yk and Y{tk). Nevertheless, Y k from <\5.8b still apparently approximates Y(tk) so that 
the overall algorithm for computing the option price (\2.12\\ -( T2.15\\ remains of weak order 
0(A + h). This can be justified by some nonrigorous arguments and this was also demon- 
strated in our numerical experiments. To obtain such a result rigorously, we need to 
conduct convergence proof without using the intermediate finite- dimensional SDEs (\3. 10§ . 



(3.24)- We do not pursue this direction in the paper. At the same time, we note that in 



all our numerical tests the scheme using A Y from (|<5. 7p gave more accurate results than 



the scheme with A Y from ( 15. gp in the cases when the Ti nodes do not belong to the t-grid. 
Otherwise A Y in ( |5. 7fl and A Y in {\5.8§ obviously coincide. 



5.2 Algorithm of order 0(A 2 + h) 

In this section we use quadrature rules ( I3.9p . (13 . 191) and a short rate approximation ( 13 .211) 
of order 0(A 2 ). 

We aim at applying the standard composite trapezoid rule to the integrals Ij(s,Ti) 
in ( 13 .4p and ( 13 .7p . The trapezoid rule requires that each of the integration subintervals 
[T/ W ,s], [s,T g(s) ], [T g(s) ,T e(s)+1 ], [T^Ti] span at least two nodes on the T-grid. 
However, the integration intervals [T^, s] and [s, T g ] usually contain just a single node 
on the T-grid: T( s ) and T e(s ), respectively. We resolve this issue by applying the right 
and left rectangle rules on these two intervals, respectively. Thus, the quadrature rule 
Sj. (s, T, A) takes the form for s G [to, t*], % = i(s), . . . , T* : 

<S^(s,T?( s ), A) = (T (s) - s) cr j (s,T (s) ),(5.9) 

A i_1 

Si 3 (s,Ti,A) = (T e(s) - s)aj(s,T e ^) + — ^ [er 3 -(s,T m ) + o-j(s, T m+1 )\ for s < T, 

rn=g{s) 

j = 1, . . . , d. 

This quadrature rule satisfies the order condition ( 13. 6p with p = 2. To this end, we recall 
that left and right rectangle rules have local order two and we use them here on one or two 
integration steps only while the trapezoid rule has local order three and the composite 
trapezoid rule is of order two. 

To ensure that ( 15.31) satisfies Assumption 3.2' with q — 1, we, in particular, need 
to approximate the integral j t k+1 Sj^s,^, A)ds by Si. (t k , T; A, h) on a single step with 
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weak order 0{h 2 ). If the node T £k+1 is not between t k and t k+ i (due to (15. ip it cannot 
be more than one T-node in (t k ,t k+ i)), it is sufficient to approximate the integral by the 
left rectangle rule and put Sj.(t k , Ti, A, h) = hSj.{t k ,Ti, A), where Sj. (t k , Tj, A) is of the 
form ( 15. 91) but with a m j(t k ) instead of <jj(s, T m ). However, if T£ fc+1 > ^ then due to ( 15.91) 
we apply one integration rule on [tjt,T4 +1 ] and the other on [Te k+1 ,T g ], which causes 
loss of smoothness of the integrand S 1 /. (s, T i; A). To reach the required order 0(h 2 ), we 
construct the approximation using the following guidance: 

Si j (s,Ti,A)ds = / S / ,(s,:r i ,A)ds+ / S Ij (s,T i ,A)ds (5.10) 

(T £fc+1 - s)a 4+lJ (s) + —V£ k+uj {s) + ^ ek+1 ,j(s) ds 



k OT -k 



/tk + 1 
-h+i 
^ rtk+i l ~ l 

Jtk m= Bk+1 

(T + a T 4+i -t k + A 

-tk) 2 °4+ijW 

A ^ 

+ ^"2~ [^"m,j(^fc) + <5m+lj(£fc)] • 

m =gfc+i 



As a result §/(£&, Tj, A) in (15.31) is taken of the form: 



Sj^tfcjT^jA,/!) = hA ek+1)k a ek+1)j (t k ), (5.11) 
hA ek+1)k a ek+1 j(t k ), ifT £fc+1 < £ fc , 



A 4+i,fe - £ ^ (T 4 + ij(*fe) - A 4+i,fe+2 T^fc+xj^fc), otherwise, 

A /_ 11 

Sjj^k, A, h) = Si.(t k} T Qk+l ;A, h) + h — cr ek+1 ,j(tk) + 2 cr m ,j{tk) + &i,j(t k ) 

\ m=e k+1 +l 

i = Q k+ i + l,---,N, j = l,...,d. 

By a similar reasoning used to derive ( 15. 9p . we obtain the corresponding quadrature 
rule Sz(t*,T*, A) (see (13 . 19[) ) . Namely, we apply the right-rectangle rule on the integration 
interval [t M , T e ] and the composite trapezoid rule on the rest of the integration interval, 
i.e., 

s z {t\T\A) = firA eMM + ^[fir + 2 Ht+fS)- (5-12) 

It is not difficult to show that the combination of rectangle and trapezoid rules used for 
deriving (I5.12p satisfies the order condition (I3.20p with p = 2. 
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We use linear interpolation for the short rate in (I3.2ip : 
i(t*) r 



7T 



(*>=£ 



1=0 



t -^f^T l+l ) + 7 ^-^f(t } T l ) 



Xt£[T h T l+1 )i 1 e [*0,** 



(5.13) 



The approximation ( 15 .131) obviously satisfies the order condition (13.221) with p = 2. As 
in the case of Algorithm 5.1, the coefficient in the right-hand side of (13.241) is also only 
piece-wise smooth here. Consider first the case when the node is not between 

t k and t k+1 . The application of the left-rectangle rule to the integral f* k+1 n(s)ds has 
the error 0(h 2 /A), i.e., it does not lead to a uniform error estimate 0(h 2 ) required by 
Assumption 3.2'. To ensure that the estimates 0(h 2 ) in Assumptions 3.2' are uniform in 
A, we use the following guidance: 



7i(s)ds 



tk 



rp rp 

s ~ fc+1 ./^+i( S ) + . 8k+1 



f 0k+1 (t k 
f gk+1 (t k )h 



A 

s-T, 



A 



ds 



k + l 



A 



•fc+i T — s 

1 Ck+1 b 



A 



-ds 



tk+i/2 - T t k+± + ^ k+1 ^ h ^ek+i ~ *fc+i/a 



So, in this case we put A Y (t k ;h) = h 



A 



V fc + 1 ,fc+l/2 fik+1 



A 



ft-k+1 _£fc+l»fe+Vg f Sk+1 

Jk A -Ik 



In the other 

tk 

k + l 



case, i.e., if Tg k+1 > t k) we split the integral f t * +1 n(s)ds = J tk k+1 n(s)ds + 7r(s)ds 

and approximate each of them separately as we did in constructing (15. 7p . As a result, we 
arrive at 



A Y (t k ;h) = { 



h 
A 



A P k+1 ,k+l/2 jl k +l 



■ ft-k + l 
Jk 



V fc+1 ,fc+l/2 rSk+ 
Jk 



•fc+li 



A J k 
2 A J k 2A J I- 



A e k +i^ Ilk 
Jk 



if T t ... < t k , 



Jk 



(5.14) 



-A 



e k +i,k+i 



A g fc+1 +i.fc+i itk+i _ A ^fc+i,fc+i IQk+i 
2A J k+l 2A Jk+1 



otherwise. 



Assumption 3.2' with q = 1 can be checked for the scheme (15. 3ft . (15.111) . ( 15.14ft following 
the standard way. 

The algorithm based on (1573]) and ([5TTTj) . (1512]) . (I5TT4|) we will call Algorithm 5.2 
for the option price (I2.12p - (l2.15p . According to Theorem I4.4[ this algorithm is of order 
0(A 2 + h). In practice (see Remark 14. 3[) we choose A = ay/h with a > such that 
(15. ip is satisfied, which results in the algorithm's accuracy 0(h). In our experiments (see 
Section [7]) Algorithm 5.2 outperformed Algorithm 5.1. 



5.3 Algorithm of order 0(A 4 + h) 

At the beginning of Section [3711 we made the assumption that there is a sufficient number 
of nodes Tj between t* and T* which ensures that we have enough nodes on the T-grid for 
using the quadrature rules ( 13. 9 p and (13 . 1 9[) and the short rate approximations (I3.2ip of 
the required accuracy. This assumption gives an unnecessary restriction for using higher- 
order algorithms in practice and we now demonstrate how it can be relaxed. To this end, 
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we introduce N instead of N in the method (13. 1 5 p as the number of discretization nodes 
on T-grid: 

N' := N V max K(f,TA V (£(t*) + 9) , (5.15) 

0<i<N 

where «(t*,Tj) and 9 are as in (13. 9p and (I3.2ip . respectively. Also, in (13 . 19[) we can put 
N instead of N and if required increase N further to be able to approximate the integral 
Z(t*,T*) on the left-hand side of (13.191) with the prescribed accuracy. As a result, we 
avoid the restriction on how close t* can be to T*. It is clear that this extension of the 
T-grid by a fixed number of nodes in the case of large A does not influence our theoretical 
results. 

Without re- writing the Euler-type scheme ( 15.31) . we will assume in this section that 
we run it for i = £k+i, . . . ,N instead of i = £k+i, ■ ■ ■ ,N. 

We are aiming at constructing an algorithm of order 0(A 4 + h) and would like to 
exploit the standard composite Simpson rule for approximation of the integrals Ij{s, T) = 
f 1 <7j(s,u)du from (13.41) and (13.71) . The Simpson rule needs three nodes per integration 
step. But the integrals Ij(s,Ti^), Ij(s,T g ( s )), and Ij(s, T e ( s )+i) are over the intervals 
which have just one or two nodes on the T-grid under (15.11) . 

We first consider the integrals Jj(s,T) = fj'<jj(s,u)du with T = T^ s ), T Q ^, and 
Tq( s )+i, which we approximate by quadrature rules S/.,(s,T, A) of the form 

S^Ti, A) = (T - s) [Pivj&Ttu) + /3* 2 <7 J -(s,T cW ) +^(s,T ff(8)+1 )] , (5.16) 

where the coefficients (3\, /3\, fl\ depend on the value of Tj. We require that (15.161) is of 
order 4, i.e., that (I3.6P is satisfied for these three integrals with p — 4. One can show that 
the following sets of coefficients satisfy this order requirement: 

Ms) _ 5 5 T e{B) -s 1 (T e(a) -sf 

Pl ~ 12 + 12 A + 6 A 2 ' {b - U) 

e(s) = 2 _ 1 T e(a) - s _ 1 (T g(s) - s) l{s) = 1 1 Tgjs) ~ s 1 (T g(s) - s) _ 

3 3 A 3 A 2 ' 3 12 12 A 6 A 2 

,2 .. /„ \2 



R g(s) _ 1 T e{a) ~ s 1 ( T 9{s) ~ s) e{s) _ 1 (T g{s) - S y 
Pi ~ i^A^ + e A 2 _i 3 A 2 ' {bAt 



Ps 4 A 6 A 2 



R s(s)+i _ 1 1 T g(a) - s 1 (T g(a) - s) J 

^ -"12 + 12 + 6 A 2 ' (5 ' 19) 



g(s)+ i = 2 1 T e(8) - s _ 1 (T e(8) - s) 2 _ (s)+1 = _5_ _ 5 r o(8) - s 1 (T e(8) - s) 2 
P * 3 3 A 3 A 2 ' Ps 12 12 A 6 A 2 

Further, for g(s) + 1 < % < N we write Jj(s,T) = Ij(s,T g ^) + Ij(T e ( a ),Ti) s) with 
Ij(T g ( s ),Ti; s) := <jj(s,u)du; and we approximate the integral Ij(T g ^, T; s), £>(s) + 
1 < i < iV , by the composite Simpson rule Sj^Tg^, T i: A; s) if its integration interval 
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spans an odd number of maturity time nodes: 

A / (i- e {s))/2-\ 
Si s ^T {,),Ti, A;s) = — [a j (s,T g{s) ) + 2 £ ^(s, T ffW+al ) (5.20) 

■j y z=i 

(*-eW)/a 

+4 XI o-j(s,T e(s)+2 /_i) + cr^s,^; 
z=i 

and otherwise we apply the Simpson's 3/8 rule for the last four nodes: 

A / (i-ff(«)-i)/2-a 
5/.(T e(fl) ,Ti,A;s) = — £Tj(s,T e ( s) ) + 2 £ ^-(s, T e(s)+2 z) 

<j y «=i 

(i-e(s)-i)/2-i 

+4 X < 7 i (s,T c(a)+2i _ 1 ) + < 7 i (s,T^ 3 ) ) (5.21) 
z=i 

3A 

+ (°"i( S ) T i-3) + 3<Tj(s, T;_ 2 ) + 3c7-j(s,Ti_i) + <Tj(s,Ti)) . 
o 

By straightforward calculations one can show that the quadrature rule ( 15. 16ft . (15.201) . and 
(15.2 ip satisfies the order condition (13. 6 j) with p = 4. 

To obtain Sj. (4, T; A, zz) based on (I5.16p . ( I5.20p . and (I5.2ip . we need again to consider 
the two cases: when 4+i = 4 and hence T k+1 < tf- and when (see also ( 15.11) ) 4+1 = 4 + 1 
and hence T £k+1 > t k . If 4+1 = 4, then we put 8^.(4 ,^; A, h) := hSi.(t k ,T h A;t k ) with 
SiAtj-jTi, A;t k ) having the form (15.161) . ( I5.20I) - (I5.21I) with cri 7 j(t k ) instead of (jj(tk,Ti). 
Otherwise, we split the integral 

Sj^s^Ti, A)ds = / ^.(s,^, A)ds+ / (s, T h A)ds 

Jt h JT lk+1 

and approximate each of them to obtain the required Sj. A, h) analogously to how 

we have proceeded in constructing Algorithm 5.2. We do not write the corresponding 
expressions of §j. (t k , T; A, h) here though there is no difficulty to restore them. 

Using (I5.16p . (I5.20p . and (15.211) . we construct the quadrature rule Sz(t* ,T* , A) (see 
(I3.19P ) and arrive at 

S z (t\ T, A) = A iM [ftfft + tifl? + P 3 flr +1 ] ; (5-22) 

a / (n- 8m -i)/2 _ (N-g M +i)/2 

Sz(T 8M ,T\ A) = - f*? + 2 £ f^ +21 + 4 £ /J^ 1 + ] ; (5.23) 

<j y 1=1 1=1 

A f - ( Ar -<?Af)/2-2 _ [N-g M )/2-\ 

S z (T eM ,T*,A) = - /If +2 E /lf +2 ' + 4 £ /lr + +S" S 

y /=i z=i 

+ ^ 3 + 3/r 2 + 3/Hf- 1 + ffi) . (5.24) 
Then we define Sz(t* , T*,A) to be used in the algorithm as 

(IE22D, dSHBD with i = qm if iv= 

S z (t*,T*, A) = <( (JE22]), d539D with i = £ M + 1 if iV = £ A/ + 1, (5.25) 
S z (t M , Tg u , A) + S z (T eM ,T N , A) if N > g M + 1, 
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where Sz(t M ,T 0M , A) is from (15.221) . (I5.19P with i = g M and 

a Z [i QM ,i N ,is) - i (JS23 j ) ifN _ QM + liseyen 

It is not difficult to show that the quadrature rules used for deriving f !5.25|) satisfy the 
order condition (I3.20p with p = 4. 

For the short rate approximation rr(t) (see (I3.2ip ). we use cubic polynomial interpo- 
lation which obviously satisfies the order condition (I3.22p with p = 4: 

3 

7r(t) = J2L 3 (t)f(t,T m+J ), (5.26) 

j=0 



where 

L&) = n 



t — T£(t)+i 



To obtain the corresponding A (tk',h) = A fl,j = Ik, ■ ■ ■ , ^k+i +3; ft,), we follow 
similar guidance as the one used to obtain (15.141) . We do not write the expression of 
A Y (tk] h) here but there is no difficulty to restore it. 

The algorithm presented in this section, we will call Algorithm 5.3 for the option 
price (I2.12p -f l2.15p . Assumption 3.2' with q — 1 can be checked for this algorithm following 
the standard way. According to Theorem 14.41 Algorithm 5.3 is of order 0(A 4 + h). In 
practice (see Remark l4.3p we will choose A = a\fh with a > such that (15. ip holds, 
which results in the algorithm's accuracy 0(h). 

Remark 5.2 (Complexity of the algorithms) Let us estimate computational complexity of 
the algorithms considered in this section. The number of operations in these algorithms is 
of order O(MN). Then running times of Algorithm 5.1 with A = ah, Algorithm 5.2 with 
A = a\/h, and Algorithm 5.3 with A = a\fh are proportional to M 2 , Ma/M, and M\[M, 
respectively. Also, it should be taken into account that Algorithms 5.2 and 5.3 require 
approximately twice and four times number of operations per t-step, respectively, than 
Algorithm 5.1. Hence one can expect that in reaching a similar accuracy Algorithm 5.3 is 
approximately M 3 / 4 /4 faster than Algorithm 5.1 and Algorithm 5.2 is yM/2 faster than 
Algorithm 5.1. This is confirmed in our numerical experiments (see Section^ . 

Remark 5.3 // in Algorithms 5.2 and 5.3 we substitute the Euler scheme (\5.3\i by a 
second-order (i.e., q = 2) weak scheme (see examples of such schemes in, e.g. \TBj ), then 
these modified algorithms (they should satisfy Assumption 3.2' with q = 2) will become of 
order 0(h 2 + A 2 ) and 0(h 2 + A 4 ), respectively. Choosing A = h and A = ay/h in the 
modified Algorithms 5.2 and 5.3, respectively, their accuracy becomes of order 0(h 2 ). 



6 Mean-square method and its convergence 

In most of the financial applications weak numerical methods, which we have considered in 
the previous sections, are sufficient. At the same time, mean-square methods can be useful 
for simulating scenarios. Also, mean-square convergence of fully discrete approximations 
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for the HJM model is of theoretical interest. In this section we consider a mean-square 
method for (I3.10p - fl3.12p and prove its convergence. 

We consider an approximation fl +1 of f l {t k+ i) from f)3.10p (i.e., a full discretization 
of f|2T4jl - (j23|l in both T and t) of the form 

f = f , i = 0,...,JV, (6.1) 

fk+l = /fc+ 

+A%, T<; j = 4+i, • • • , K(t fc+ i, Tj) V z; /i; - W t (t k ), I = 1, . . . , d, t k < s < t k+l ), 

i = £ k+1 ,...,N, k = 0,...,M, 

where the form of function A 1 depends on the coefficients of fl3.10D - fl3.12p . i.e., on a and 
on choice of the quadrature rule SY,; K-(t k ,Ti) is as in the quadrature f)3.12p . Note that in 
this section we use the same notation f l k for the mean-square approximation as the one 
we use for weak approximations in all the other sections of this paper. Since mean-square 
approximations of (13.101) are considered in this section only, this abuse of notation does 
not lead to any confusion. 
As before, we put 

n = fL, fc = m + l,...,M, 0<i<£{t*)-l, 

where m = [(T i+1 — t ) /h~\ — 1. Then the N + 1-dimensional vector {fl, i = 0, . . . , iV} is 
defined for all k = 0, . . . , M. 

We impose the following assumption on the one-step approximation fl x (t + h) of the 

method (16. ip for the solution fl x {t + h) of f )3.10p with the initial condition x given at time 

t:fl x (t)=x\ 

Assumption 6.1 Let 

Q2 > ~ , qi > Q2 + ~ ■ (6.2) 

Suppose the one-step approximation ]\ x {t + h) has order of accuracy q± for expectation 
of the deviation and order of accuracy q<i for the mean-square deviation; more precisely, 
for arbitrary t < t < t* — h, x G IR 7V+1 the following inequalities hold: 

\E(fUt + h)-fUt + h))\<Ch^, (6.3) 

E\fUt + h)- fl x (t + h))\ 2 ] V2 < Ch«\ (6.4) 
i = 0,...JV, 

where C > is a constant independent of h, A, and x. 

Assumption 6.1 is analogous to the conditions of the fundamental theorem of mean- 
square convergence [TH p. 4]. We note that C in fl6.3p - fl6.4l) are independent of x while 
in the fundamental theorem such C depend on x. In our case it is natural to put C inde- 
pendent of x since the coefficients of f)3.10p and their derivatives are uniformly bounded 
(see Assumptions 2.1-2.2). We also emphasize that the constants C in fl6.3p - fl6.4D do not 
depend on A. 

Under the stated assumptions we will prove mean-square convergence of f k to f l {t k ) 
uniform in A in order then to prove mean-square convergence of f % k to f(t k , Tj) exploiting 
in addition Theorem 14. 11 We cannot use here the fundamental theorem of mean-square 
convergence [T5| p. 4] since we need to show that the convergence is uniform in A. 
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Theorem 6.1 Suppose Assumptions 2.1-2.3 and Assumption 6.1 are satisfied. Then for 
any M, N and k = 0, 1, . . . , M, i = 0, 1, . . . , N the following inequality holds: 

- "I 1/2 

E\f(t k )-fi\ 2 \ <Kh^ l '\ (6.5) 
i.e., the order of mean-square accuracy of the method (1 6. II) for ( \3. l(k is q = q2 — 1/2. 
Proof. We have (cf . (HJ pp. 7-8]) 

f(t k+1 ) - fl +1 = ft ,f (tk+i) - ftojoitk+i) = ft k J(t k )( tk +^ ~ ft k jS tk+1 ^ ( 6 - 6 ) 

= ^Lmfik+i) - flj^k+i)) + (/t k ,/ k (*fc+i) - fLiS^+i)) , 

where the first difference in the right-hand side of (I6.6P is the error of the solution arising 
due to the error in the initial data at time tk, accumulated over k steps, and the second 
difference is the one-step error at the (k + l)-step. Taking the square of both sides of 
(16.6p . we obtain 

Rl k+1 :=E\f(t k+1 )-ft +1 \ 2 (6.7) 
= ^(l^)^)-^(Wi)| 2 |^) 

+^(l^(**+0-^(**+i)l 2 I^O 

Due to the condition ( 16. 4p . we get for the second term on the right-hand side of ( 16.7[) : 

\EE(\fl Jk (t k+l ) - fl Jk (t k+l )\ 2 \T tk )\ < Ch 2 <*. (6.8) 

Let us estimate the first term on the right-hand side of ( 16. 7p . Ito's formula implies 
that 



£ i(tk+l) ■- E \ft k ,f(t k )( tk +^ ~ ft h ,fS tk+1 



ci 1 2 



+2£ /' +1 ^j(..»W-4,/.W) 

Jt k 

x [a T (s, T h £ kj/(tfc) f(t k ); s, T h A) - a T (s, T h f\ kh {s))~S l {t k , f k ; s, T u A)]ds 

+E \a{s,T h fl kJ{tk) {s)) -a{s,T h fl kJk {s))\ 2 ds. 

J th 



Then, recalling that a(s,T,z) is globally Lipschitz in z due to Assumption 2.2 and the 
form of Si(t k , f k ; s,Ti,A) (see ( I4.35P ). we obtain 



eKh + i) < E\f\t k ) - fl\ 2 + K r +1 ^l4j ( , fc) (^)-i5 fc ,/ fc (^)| 2 ^ 



l=£{s) 
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where K > does not depend on A. Introduce £ 2 Max (s) := max 0< . <M s 2 (s). Then 



4«M < maX - fk\ + K / 4/ax(^)^ 



; fc+l 



which implies that for all < i < M and all sufficiently small h > : 

e?(f*+i) < ^ max £|f (t fe ) - /^| 2 < max £|f (t fe ) - /*| 2 • (1 + Kh), (6.9) 



where > does not depend on A and h. 

Now let us estimate the third term on the right-hand side of (16.71) . We have (cf. 
Lemma 1.1.3 in [T8l p. 5]): 

fim^+i) ~ fU^+i) = Hh) -ft + z\ (6.io) 

where 

"*fc+i 



jf W [a T (s, Ti, j? fc)/(tt) (*))&(**, f(t k ); s, T it A) 
-a T (s, T i? (s))5r(t fc , /*; s, T*, A)]ds 



Using (16. 9p . it is not difficult to get 

E (Z*) 2 < Kh ■ max E\f(t k ) - f k \\ (6.11) 

0<i<M 

where K > does not depend on A and Using flBTTO]) . (Oil . (IO> . (EHIjl . and 
we obtain 

- - ( 6 - 12 ) 

< |£(r(^) - f k )E{fi kh {t k+l ) - r tkJk (t k+1 ))\F th )\ 

< (^if(4) - m 2 ) 1 ' 2 ■ Kh* + (e {zf) 1/2 (E(fi Jk (t k+1 ) - r tkJk {tk + if) 1/2 

( _ \l/2 

< Kh*(E\f(t k ) - ill 2 ) 1 / 2 + Kh^' 2 ■ max E\f(t k ) - f k \ 2 ) 

\0<i<M J 

/ \ V2 

< ^ +1 / 2 - max£|f(4)-/*| 2 , 

\0<i<Af J 

where K > does not depend on A and h. 

Let i^ aXifc := max Q <.< M # 2 fc . Then it follows from (iO) . flETSft . (P and (15321) that 

< i^aM • (1 + + Kh^ 2 R Max , k + Ch 2 ^. 
Using the elementary relation 



h^ 2 R Max , k < + h 
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we get 

^Max,k+l — RMax,k ' (1 + ^h) + Ch 2q2 

whence (16.5p follows taking into account Lemma 1.1.6 from [18, p. 7] and the fact that 
R 2 Max0 = 0. Theorem 16. II is proved. □ 

Theorems 14.11 and 16.11 imply the following result. 

Theorem 6.2 Assume that the conditions of Theorems 4-1 and \6.1\ hold. Then for any 
M, N and i = 0, 1, . . . , N, k = 0, 1, . . . , |~(T, — t ) /h] — 1 the mean-square error is esti- 
mated as 

[E\f(t h ,T t ) - ft\ 2 ] 1/2 < K ■ (A p + h^ 1 ' 2 ) , (6.13) 
where K > is a constant independent of A and h. 

Example 6.1 To illustrate the results of this section, let us present a mean-square algo- 
rithm for ( |#.^p - (lff.5p based on the mean-square Euler-type scheme: 

r = r , i = o,...,N, f = o, (6.i4) 
fi+i = f k + e m**)Sj, a, h) + h 1 ' 2 £ ^ij(**)e ilfc+ i, 

i=i i=i 
i = 4+i) • • • ,N, 

where £ • fc+1 are independent Gaussian random variables with zero mean and unit variance, 
(vi,i(t k ), ■ ■ ■ , o"i,d(4)) T = (ax(tk,Ti, fl), . . . ,a d (t k ,Ti, f k )) J , 

and Sj. (4> Ti, A, h) depends on our choice of the quadrature rule (jff.ffp . If Sr. (4, T»; A, /i) 
zs taken from { 5.4} or (15. 1 il) or from Algorithm 5.3 then p = 1, p = 2 orp = 4, respectively, 
and qi = 2 and q% = \. under h < A. The overall error of these algorithms are 0(A + h 1 ^ 2 ), 
0{A 2 + h 1 / 2 ), and 0(A 4 + h 1 / 2 ), respectively. 



7 Numerical examples 

In this section we demonstrate accuracy and convergence properties of the algorithms 
from Section O We also compare computational costs of the algorithms. This comparison 
illustrates that the algorithms with higher-order quadrature rules are more efficient. 

For illustration, we price an interest rate caplet which is an interest rate derivative 
providing protection against an increase in an interest rate for a single period. Suppose 
a caplet is set at time t* with payment date at T* and has the unit nominal value and a 
strike K. The arbitrage price of the caplet is given by (12. 9p with t* = s k and T* = Sj. The 
caplet parameters chosen for the experiments are t* = 1.0, T* = 6.0, K = 0.03. 

A particular model within the HJM framework (I2.4p - (I2.5I) is specified by a choice of 
the volatility function and initial forward rate curve. Here we consider two examples: a 
one-factor model with deterministic exponential volatility function (Vasicek model, see, 
e.g. H]) and a two-factor model with proportional volatility function (see, e.g. [131 El 
[2T],[lTJ). The former one admits a closed- form formula for the caplet price. 

The algorithms were implemented using C++ with GCC 3.4.3 compiler. The exper- 
iments were run on ALICE HPC Computer nodes of the University of Leicester, each 
with dual quad-core 2.67GHz Intel Xeon X5550 processor, 12 GB RAM, and OS 64-bit 
Scientific Linux 5.4. 
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7.1 Vasicek model 



We consider the one-factor HJM model fl2.4p - fl2.5l) with the deterministic volatility func- 
tion given by 

o~(t, T) = cr exp( — k(T — £)), (7.1) 
and the initial forward curve defined as 

a 2 , 

fo T) = exp(- K (T - t r + (1 - exp (-k T - t ))) - r-^ 1 - exp (-« T - *„ , 

(7.2) 

where a, k, r , and $ are given positive constants. 

Table 7.1: Algorithm 5.1 for the Vasicek model. Performance of Algorithm 5.1 with 
A = h in the case of the Vasicek model (17. ip . (17.21) with parameters a = 0.02, r = 0.05, 
k = 1 and 9 = 1 for pricing a unit nominal caplet with parameters to — 0, t* = 1.0, 
T* = 6.0, K = 0.03. L is the number of independent runs in the Monte Carlo simulation 
(see (ESP). 



h 


L 




error 




CPU time, min 


0.2 


10 7 


4.22 x 10" 


- 2 ±2.80 x 10~ 


-6 


4.00 x 10" 1 


0.1 


10 7 


2.04 x 10" 


- 2 ±3.06 x 10" 


-6 


9.00 x 10- 1 


0.05 


10 7 


1.00 x 10" 


- 2 ±3.19 x 10- 


-6 


2.45 x 10° 


0.025 


10 7 


4.98 x 10" 


- 3 ±3.26 x 10- 


-6 


7.73 x 10° 


0.0125 


10 9 


2.48 x 10" 


- 3 ±3.29 x 10- 


-7 


2.30 x 10 3 


0.00625 


10 9 


1.24 x 10" 


" 3 ± 3.31 x 10" 


-7 


8.32 x 10 3 



Table 7.2: Algorithm 5.2 for the Vasicek model. Performance of Algorithm 5.2 with 
A = y/h in the case of the Vasicek model ( 17. ip . (17. 2p with the same parameters as in 
Table EU 



h 


L 




error 




CPU time, min 


0.2 


10 7 


6.53 x 10' 


- 3 ± 2.72 x 10" 


-6 


3.00 x 10- 1 


0.1 


10 7 


3.32 x 10" 


- 3 ±2.99 x 10- 


-6 


5.33 x lO" 1 


0.05 


10 7 


1.65 x 10" 


- 3 ±3.11 x 10- 


-6 


1.03 x 10° 


0.025 


10 7 


8.29 x 10" 


- 4 ±3.22 x 10" 


-6 


2.13 x 10° 


0.0125 


10 9 


4.13 x 10" 


- 4 ±3.27 x 10- 


-7 


4.65 x 10 2 


0.00625 


10 9 


2.07 x 10" 


' 4 ±3.30 x 10- 


-7 


1.09 x 10 3 



It is known (see, e.g. [2j HJ [23]) that a caplet corresponds to a put option on a zero- 
coupon bond. In [TJ] analytic expressions for the European option prices on zero-coupon 
and coupon bearing bonds under the Vasicek model are derived. In particular, the price 
of the caplet set at time t* with payment date at T*, unit nominal value and strike K is 
given by 

F(to, f (■) ; t*, T*) = P(t , t* )$(-c P + ap) - (1 + K(T* - t*))P(t , T*)$(-c P ), (7.3) 
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Table 7.3: Algorithm 5.3 for the Vasicek model. Performance of Algorithm 5.3 with 
A = a\fh in the case of the Vasicek model (17. ip . ( 17. 2p with the same parameters as in 



Table ED 



h 


L 




error 


CPU time, min 


a 




0.2 


10 7 


1.25 x 10" 


- 3 ±2.53 x 10~ 6 


2.59 x 10- 1 


9.97 x 10- 


-i 


0.1 


10 7 


6.28 x 10" 


" 4 ±2.87 x 10" 6 


4.10 x lO" 1 


9.70 x 10" 


-i 


0.05 


10 7 


3.18 x 10" 


" 4 ± 3.09 x 10" 6 


8.11 x 10" 1 


9.76 x 10- 


-l 


0.025 


10 7 


1.56 x 10" 


" 4 ±3.20 x 10~ 6 


1.59 x 10° 


9.43 x 10- 


-i 


0.0125 


10 9 


9.62 x 10" 


- 5 ±3.26 x 10~ 7 


3.13 x 10 2 


9.97 x 10- 


-l 


0.00625 


10 9 


4.71 x 10" 


- 5 ±3.30 x 10~ 7 


5.96 x 10 2 


9.70 x 10- 


-i 



where $(•) denotes the standard normal cumulative distribution function and 

[l-exp(-2/c(T*-t*))], 



a 



o P 



exp(-2ft(t* - t )) 



cp 



-lln 

Op 



2k 

[1 + K(T* -t*))P(t ,T*) 

P{to,r) 



The values of parameters chosen in the experiments are t = 0, a = 0.02, r = 0.05, 
k — 1, and d — 1. The values k — 1 and d — 1 are rather unrealistic from the financial 
point of view and are chosen for illustrative purposes. Under a more realistic choice of 
parameters, simulations are done with a particular time step h (see Table 17741) . 

The results of the experiments with Algorithm 5.1 of order 0(A + h), Algorithm 5.2 of 
order 0(A 2 + h), and Algorithm 5.3 of order 0(A 4 + h) are presented in Tables 17.11 17.2[ 
and !7.3[ respectively. For Algorithms 5.1 and 5.2, we set A = h and A = V7t, respectively. 
For Algorithm 5.3, we set A = ayh with a > so that T-grid remains equally spaced: 



a 



( T*-t 


1 


T*~t 






1 T*-t 




T*-t Q 




^ J 



+1 



otherwise, 



T*-t 



> 



i 

2' 



(7.4) 



where |_ - J denotes the integer part of a real number. It is clear that a ~ 1. 

As a result, the errors of all three algorithms become of order O(h). In the tables, 
the values before "±" are estimates of the bias computed as the difference between the 
exact caplet price ( 17. 3p and its sampled approximation (see ( 13.34(1 ). while the values 
after "±" give half of the size of the confidence interval for the corresponding estimator 
with probability 0.95. The number of Monte Carlo runs L is chosen here so that the 
Monte Carlo error is small in comparison with the bias. It is not difficult to see from 
the tables that the experimentally observed convergence rate is in agreement with the 
theoretical first order convergence in h. We note that in the analysis of convergence of 
Algorithm 5.3 one has to take into account not only values of h but also of a. As expected, 
the experiments demonstrate that Algorithm 5.3 is the most computationally efficient 
among the three algorithms tested and also Algorithm 5.2 outperforms Algorithm 5.1. As 
it follows from Tables I7.1ti7.4[ for a fixed time step h the ratios of running times of the 
considered algorithms is in agreement with the theoretical prediction (see Remark 15.21) . 

In Table 17.41 we present results for h = 0.1 and L = 10 9 in the case of the more 
realistic choice of parameters k = 0.178 and d = 0.086 of the Vasicek model. With 
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Table 7.4: Vasicek model. Performance of the algorithms (I5.ip - (I5.3I) with h = 2~ 6 and 
L = 10 9 in the case of the Vasicek model (17. ip . (17.21) with parameters a = 0.02, r = 0.05, 
k = 0.178 and 9 = 0.086 for pricing a unit nominal caplet with parameters to = 0, 
t* = 1.0, T* = 6.0, K = 0.03. 





h 


L 


error 


CPU time, rain 


Algorithm 5.1 


0.1 


10 9 


-5.38 x 10" 4 ±2.90 x 10" 6 


9.79 x 10 1 


Algorithm 5.2 


0.1 


10 9 


1.75 x 10^ 4 ±2.86 x 10~ 6 


5.51 x 10 1 


Algorithm 5.3 


0.1 


10 9 


7.81 x lO" 5 ±2.89 x 10~ 6 


4.10 x 10 1 



these parameters, the bias is very small, and if one would like to analyze it, e.g. for 
h = 2~ 2 x 5 _1 , then the number of Monte Carlo runs has to be increased up to 10 11 
in order to make the Monte Carlo error sufficiently smaller than the bias. We see from 
Table EH that Algorithm 5.3 is more than twice faster than Algorithm 5.1 in producing 
the results of a similar accuracy. 

7.2 Proportional volatility model 

Here we choose the volatility functions of the form 

aj(t, T) = Gj exp(-/t i (T - t)) min (f(t, T), T) , (7.5) 

where aj and Kj are positive constants and T is a large positive number introduced to 
cap the proportional volatility in order to avoid an explosion of the forward-rate process 
(cf. Assumption 2.1 and also Remark 12. ip . The volatility specification of the form (17.51) 
yields an approximately lognormal distribution of forward rates. 

Table 7.5: Algorithm 5.1 for the Proportional volatility model. Performance of Algo- 
rithm 5.1 with A = h in the case of the proportional volatility model (17. 5p with param- 
eters (I7.6P and with initial forward curve (17. 7p for pricing a unit nominal caplet with 
parameters t = 0, t* = 1.0, T* = 6.0, K = 0.03. 



h 


L 




error 


CPU time, min 


0.2 


10 9 


5.80 x 10" 


- 4 ±2.57 x 10" 6 


6.99 x 10 1 


0.125 


10 9 


3.64 x 10" 


- 4 ±2.57 x 10" 6 


1.21 x 10 2 


0.1 


10 9 


2.92 x 10" 


- 4 ±2.57 x 10" 6 


1.67 x 10 2 


0.05 


10 9 


1.48 x 10" 


- 4 ±2.56 x 10~ 6 


4.84 x 10 2 



Let us note that in (11] a number of volatility models including one and two factors 
proportional volatility models are examined. The performance of the models is evaluated 
based on the accuracy of their out-of-sample price prediction and their ability to hedge 
caps and floors. This study reveals that in out-of-sample pricing accuracy the one- and 
two- factor proportional volatility models outperform the other competing one- and two- 
factor models, correspondingly. The one-factor BGM model outperforms the proportional 
volatility model only in pricing tests, which were not strictly out-of-sample. In terms of 
hedging performance, the two-factor models provides significantly better results than the 
one-factor models. 
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Table 7.6: Algorithm 5.2 for the Proportional volatility model. Performance of Algo- 
rithm 5.2 with A = \fh in the case of the proportional volatility model (17. 5p with the 
same parameters as in Table 17.51 



h 


L 




error 


CPU time, min 


0.2 


10 9 


1.50 x 10" 


" 4 ±2.56 x lO" 6 


5.70 x 10 1 


0.125 


10 9 


8.91 x 10" 


- 5 ±2.56 x 10~ 6 


7.88 x 10 1 


0.1 


10 9 


8.01 x 10" 


- 5 ±2.56 x 10~ 6 


1.03 x 10 2 


0.05 


10 9 


3.47 x 10" 


- 5 ±2.56 x lO" 6 


2.16 x 10 2 



Table 7.7: Algorithm 5.3 for the Proportional volatility model. Performance of Algo- 
rithm 5.3 with A = a\fh in the case of the proportional volatility model (I7.5jl with the 
same parameters as in Table 17.51 



h 


L 




error 


CPU time, min 


a 




0.2 


10 9 


7.04 x 10" 


' 5 ±2.57 x 10~ 6 


5.32 x 10 1 


9.97 x 10- 


-l 


0.125 


10 9 


4.59 x 10" 


- 5 ±2.57 x 10~ 6 


6.34 x 10 1 


9.17 x 10- 


-l 


0.1 


10 9 


3.66 x 10" 


- 5 ±2.57 x KT 6 


7.40 x 10 1 


9.70 x 10- 


-i 


0.05 


10 9 


1.74 x 10" 


- 5 ±2.57 x 10~ 6 


1.54 x 10 2 


9.76 x 10- 


-l 



In our experiments we consider two factors, i.e., d = 2. We use the same parameters 
for (I7.5P as those found in [11] by calibrating the model to the market prices of caps and 
floors across different maturities and strike rates: 

a x = 0.1043, a 2 = 0.1719, m = 0.052, k 2 = 0.035. (7.6) 

As the initial forward curve, we take the one used in numerical examples in [TO] : 

f (T) = log(150±48T)/100. (7.7) 

Since the closed-form formula for caplet price is not available for the HJM model 
(I2.4p - (l2.5p with the volatility ( 17. 51) . we found the reference caplet price by evaluating the 
price using Algorithm 5.3 with h = 0.00625, A = a\fh with a from ( I7.4p . and taking the 
number of Monte Carlo runs L = 10 9 . This reference value has the Monte Carlo error 
2.56 x 10 -6 , which gives half of the size of the confidence interval for the corresponding 
estimator with probability 0.95. 

Tables [73| 17.61 and 17.71 report the results of our experiments for Algorithm 5.1 with 
A = h, Algorithm 5.2 with A = yh, and Algorithm 5.3 with A = a\/~h, a is from 
( I7.4p . As in the previous tables, the error column values before "±" are estimates of the 
bias computed using the reference price value and the values after "±" reflect the Monte 
Carlo error with probability 0.95. As in the Vasicek model example, the Monte Carlo 
error was made relatively small in order to be able to analyze the bias. One can observe 
that the results demonstrate first order of convergence which is in agreement with our 
theoretical results. The experiments also clearly illustrate the computational superiority of 
Algorithm 5.3 whereas Algorithm 5.1 is the slowest out of the three algorithms presented. 
The computational times are consistent with the theoretical complexity of the algorithms 
described in Remark 15.21 
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